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ABSTRACT 


The  aim  of  the  project  is  to  apply  new  or  recently  developed  computer 
graphic  techniques  to  a particular  discipline  with  the  thought  ot  broad- 
ening its  research  capabilities.  The  discipline  chosen,  that  of  medicine, 

has  emphasis  placed  upon  the  area  of  hemodynamics . \ It  is  seen  that  com- 

\ 

puter  graphics  methods  lend  themselves  admirably  to /studies  of  blood  flow 

y 

subsequent  to  analyses  of  particular  problems  by  the  Application  of  fluid 

L 

mechanics  concepts  and  numerical  analysis  techniques.  ^Projects  chosen  for 
contemplation  were  the  design  of  optimized  flow  prosthetic  heart  valves 
and  the  study  of  vascular  degeneration  and  replacement.  A developed  mock 
circulation  is  to  be  employed  concurrently  for  furnishing  verification  of 


the  theoretical  resu  lts^a  relatively  uncommon  attempt  to  expe rimental ly 


verify  particular  compute*:  graphics  algorithms. 


Prosthetic  heart  valves  may  be  classified  for  fluid  flow  operational 
purposes  according  to  the  presence  of  an  unobstructed  central  fluid  path- 
way upon  valve  opening  (leaflet  type) , or  a central  obstruction  with 
lateral  flow  (ball  or  poppet  type).  In  all  situations,  blood  as  the  fluid, 
sees  obstacles  in  its  path  of  flow  and  the  project,  therefore,  is  to  attempt 
to  duplicate  mathematically  this  flow  of  fluid  about  such  an  obstacle.  In 
turn,  i f an  idealized  obstacle  shape  is  transferred  from  the  center  of  flow 
to  one  wall,  the  mathematical  problem  is  basically  the  same  but  now  an 
idealized  atherosclerotic  lesion  can  be  formed  for  study.  In  continuing 
to  show  the  versatility  of  the  project,  if  one  follows  the  blood  flow 
through  the  circulatory  system,  the  corners  of  arterial  side  channels  con- 
stitute mathematical  obstacles  since  abrupt  changes  in  the  flow  direction 
are  caused,  as  .well  as  distortion  of  the  flow  lines. 

In  attempt ing  to  apply  viscous  fluid  mechanics  theory  one  is  required 


to  employ  the  Navie r t okes  equations,  a partial  differiTifi.il  equation 


system  which  is  nonlinear,  thus  allowing  exact  solutions  in  only  a few 
special  cases.  These  equations  are  within  major  areas  of  application  of 
modern  applied  mathematics  since  they  exhibit  all  the  computational 
difficulties  of  other  partial  differential  equations  of  mathematical 
physics.  Direct  numerical  solution  of  the  Navier-Stokes  equations,  however, 
now  seems  within  reach  because  of  the  availability  of  large  memory  digital 
computers.  Subsequently,  the  application  of  physical  laws  to  hemodynamic 
studies  is  enhanced  by  the  attachment  of  a graphical  display  system  to 
the  computer. 
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p Fa-:  face 


The  authors  of  this  report  are  dedicated  to  the  task  of  effecting 
a useful  and  true  man-ma  .me  interface  by  the  expanded  utilization  of 
computer  graphics  as  a c<  -pi  I : - nt  to  the  digital  computer.  Their  own  area 
of  interest  is  centered  about  the  development  of  interactive  techniques 
for  the  solution  of  large-scale  numerical  problems.  Such  problems,  at 
the  present,  are  difficult,  if  not  impossible  to  solve  practically  by 
machine  computing  processes  and  require  the  guidance  and  direction  of  a 
skilled  human.  In  particular,  graphically-aided  techniques  for  the  prac- 
tical solution  of  nonlinear  partial  differential  equations  are  of  prime 
importance  in  this  study. 

It  is  the  contention  of  the  present  investigators  that  a large 
percentage,  perhaps  nearly  all,  of  algorithms  developed  for  computer 
graphics  usage  are  the  property  of  a relatively  small  elite  group  of 
computer  scientists.  If  so,  then  certainly  the  potentialities  of  such 
algorithms  are  lost  to  many  other  disciplines  for  application.  With 
this  in  mind,  it  was  decided  that  one  of  the  mo  e conservative  of  other 
disciplines  should  be  broached  in  an  attempt  to  make  them  aware  of  the 
innate  capabilities  of  the  computer  graphics  concept.  The  discipline 
of  medicine  was  so  chosen  since  they  are  particularly  sensitive  to  the 
inroads  of  mathematical  concepts.  The  medical  world  is  newly  aware  of 
the  contributions  of  the  computer  to  the  measurement,  monitoring  and  data 
collection  of  concern  to  the  treatment  of  critically  ill  patients,  how- 
ever, it  is  relatively  unaware  of  computer  graphics  usage  for  understand! ng 
medical  phenomenon.  To  this  end,  two  medical  investigations  have  been 
pursued;  the  design  of  an  optimized  flow  artificial  heart  valve  and  the 
understanding  of  the  formation  and  positioning  of  atherosclerotic  lesions. 


Bulk  of  the  funding  for  such  applied  studies  is  provided  by  the  National 
Institute  of  Health  under  qrant  HE  12201,  with  Advanced  Research  Project 
Agency  support  fc ’ the  interactive  techniques  on  the  area  of  basic  algo- 
rithms useful  for  application  to  other  ARPA  study  areas.  In  the  large, 
it  is  an  attempt  to  make  other  disciplines  cognizant  of  computer  aided 
design  and  graphics  usage  more  catholic. 

The  first  of  these  reports  discusses  the  numerical  and  resulting 
algorithms  required  for  calculating  the  transient  dynamics  of  a specific 
incompressible  fluid  problem.  The  approach  developed  the  two-dimensional 
results  and  allowed  formation  of  resulting  configurations  with  the  aid 
of  Univac  1108  and  PDP-8  digital  computers  plus  an  IDI  display  and 
Gerber  plotter.  The  companion  report  continues  the  mathematical  technique 
so  as  to  satisfy  free  surface  boundary  conditions  and  aid  in  three- 
dimensional  studies.  Here  the  desired  results  were  formed  with  the  aid 
of  PDP-9  and  PDP-10  assemblies,  plus  a Univac  1559  display  system.  The 
project  is  a continuing  one. 
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CHAPTER  I 


INTRODUCTION 

The  aim  of  the  project  that  is  presented  in  this  report  is  to  cnlargt 
rhe  usaqe  of  present  computer  qruphic  capabilities  in  a new  area  of  endoavtu 
name ly , medicine.  The  employment  of  fluid  mechanics  concepts  is  onurrent. 
In  attempting  to  apply  viscous  fluid  mechanics  theory  one  is  required  to 
study  the  Navier-Stokes  equations,  a partial  differential  equation  system 
which  is  nonlinear,  thus  allowing  exact  solutions  in  only  a few  special 
cases.  These  equations  are  within  major  areas  of  application  of  modern 
applied  mathematics  since  they  exhibit  all  the  computational  difficulties 
of  other  partial  differential  equations  of  mathematical  physis.  Algorithm 
and  concepts  formed  within  this  investigation  should  then  b<  capable  ,f 
appl 1 cabi 1 1 ty  in  other  investigative  areas. 

Although  only  specialized  situations  have  been  fully  explored  i ..  r v 
past,  of  concern  to  the  Navier  Stokes  equations,  it  appears  that  dire  t 
numerical  solution  of  such  a system  now  is  within  reach  because  of  the  avail 
ibility  of  large  memory  digital  computers.  Subsequently,  the  application 
of  physical  laws  to  hemodynamic  studies  is  enhanced  by  the  attachment  of  a 
graphical  display  system  to  the  computer. 

Two  particular  applications  are  presented  in  this  study.  )ne  i.  an 
attempt  to  study  blood  flow  through  artificial  heart  valves  and  employ  :;u  h 
information  to  eventually  design  an  optimized  flow  prosthetic  heart  valv  . 

Tin  other  application  is  the  study  of  vascular  degeneration  r tu>-  itti 
to  amplify  knowledge  of  concern  to  atherosclerosis.  Background  knowl-  Ig^  , 
of  concern  to  both  topics,  are  discussed  somewhat  in  depth,  fui  in  i|  ; • ■ . 1 1 1 • i ; 
to  interdisciplinary  knowledge  dispersion,  such  an  action  is  thought  nec- 


essary. 


tpp]  iti  n ippear,  .it.  first  jlan  e , t<  be  {uite  dissimilar. 


ver,  : it  srt  ficial  heart  valves,  blood  as  a fluid,  se<  . bst  i l< 
flow  | . i procedure,  therei  ire,  is  an  at  tempt  to  mathe- 

iupli  ite  the  flow  of  fluid  about  such  an  obstacle.  In  turn, 
if  in  1 lulized  obstacle  shape  is  transferred  from  the  center  of  flow  to 
■n*  wall,  the  mathematical  problem  is  basically  the  same  but  now  an 

■ ithi  ■ Lerot ic  lesion  is  formed  for  study.  In  continuing  to 
. a t versatility  t the  approach,  if  one  follows  the  blood  flow  through 
• Li  ulatory  system,  the  corners  of  arterial  side  channels  constitute 
;•  nemuticul  obstacles  since  abrupt  changes  in  the  flow  direction  are 
sused,  i-  well  is  dist  rti  in  of  the  flow  lines. 


\ 


i 


i 


Cl  I APT! ;r  II 


HEART  VAI.VE  MECHANICS 

The  human  heart,  for  all  its  fascinating  capabilities,  is  basically 
just  a four  chambered  pump.  It  is  located  in  the  lower  central  portion 
of  the  chest  cavity,  just  behind  the  breast  bone  and  between  the  right 
and  left  lungs.  The  male  human  heart  is  approximately  four  and  one  half 
inches  long,  three  and  one  half  inches  wide  at  its  widest  point,  and  two 
and  one  half  inches  thick.  The  male  human  heart  weighs  from  10  to  12 
ounces  while  the  female  heart  weighs  from  8 to  10  ounces.  There  is  a 
partition  that  divides  the  heart  into  right  half  and  left  half  sections. 

In  turn,  each  half  is  divided  into  two  more  chambers  which  are  termed  the 
atrium  and  ventricle.  The  atria  act  as  collecting  chambers  for  blood 
returning  to  the  heart  from  the  lungs  and  the  veins  of  the  body.  The 
atria  fill  the  ventricles  with  blood  before  each  beat  when  the  ventricles 
are  relaxed.  '1  he  ventricles  are  muscular  chambers  which  initiate  the 
circulation  of  blood  through  portions  of  the  body  and  ultimately  back  to 
the  heart.  The  left  ventricle  is  the  true  pumping  muscle  of  the  heart  since 
it  not  only  sends  blood  through  the  body  but  it  must  pump  against  gravity 
so  as  to  forward  blood  to  the  head  and  other  body  sections  above  the  heart 
region.  The  right  ventricle  needs  only  to  furnish  blood  to  the  lungs. 

Each  ventricle  has  an  inflow  and  outflow  valve.  As  shown  in  figure  1,  the 
left  ventricle's  inflow  valve  is  the  mitral  valve  and  the  aortic  valve 
is  the  outflow  valve.  The  inlet  valve  for  the  right  ventricle  ls  the  t i i — 
cuspid  valve  while  the  outflow  valve  is  the  pulmonary  valve.  A complete 
discussion  of  the  anatomy  involved  is  not  in  order,  however,  the  above 
short  discussion  introduces  the  reader  to  a brief  description  of  the  cham- 


bers that  require  valves  for  proper  blood  motion. 


4 


A-.fir  PULMONARY 
Volve  VALVE 

PULMONARY  ARTERY  \ 

to  Right  lung \ ' ^ 

SUPERIOR  — 


. aorta 


VENA  CAVA 

PULMONARY  VEINS 
TO  RIGHT  LUNG C 


TRICUSPID  VALVE 


RIGHT 

VENTRiCLE  —\ 


- PULMONARY  ARTFRY 
TO  LEFT  LUNG 


PULMONARY  VEINS 
FROM  LEFT  LUNG 


Right  W-j/Su  „:r ; '•('  X /\  ~~  Ml  I RAL  valv 

u,'um  ) [ H\\  vrvZ>-\-*\ — V 

^sfynd  V«  4, 


\\ — CHORDAE 

TENDINEAE 
LEFT  VENTRICLE 


/ 


INFERIOR  VENA  CAVA 


PAPiLLARY 

MUSCLE 


;TT  N\  v'  >\\  -yA  r i/j 

\K  \ ) ) y / VENTRICULAR 

V/  SEPTUM 


HEART  MUSCLE 


Fig.  1:  The  Chambers  of  tho 

Human  Heart 
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In  turn,  to  understand  subsequent  ~luid  dynamic  and  computer 
iraphit  requirements,  tlie  four  heart  valves  will  bo  discussed  in  gross 
detail.  Between  each  atrium  and  ventricle  there  is  an  opening  encap- 
uldte  1 by  a fibrous  ring.  Tho  valves  are  fibrous  tissue  in  leaflet 
form  which  extend  from  the  ring  into  the  ventricle.  Delicate  cords  run 
from  the  free  edges  of  the  leaflets  to  certain  muscles  of  the  ventricles 
These  rords  form  a restraining  mechanism  so  that  when  the  ventricle  col- 
1 ips"s,  these-  leaflets  close  tight  but  do  be  mi  somewhat  in  a reverse  fas 
thi  it  rium  and  allow  so  mi  L<  ikagi  . The  miti  * : valve , be  tw<  sen  the 
left  itrium  anti  ventricle  has  two  leaflets  that  act  as  a one  way  valve. 


ill  ion 


.1.'  t w ■ leaflet  ; are  uneven  in  size  and  one  >verlaps  the  other  slightly. 


In  the  case  of  the  tricuspid  valve,  which  is  the  one  way  valve  bo two  i. 
the  right  atrium  and  right  ventricle,  it  consists  of  three  leaflets  or 
cusps.  The  pulmonary  valve  allows  one  way  blood  flow  from  tit  right 
ventricle  to  the  pulmonary  artorv  while  the  aortic  valve  allows  one 
way  blood  flow  from  the  left  ventiicle  to  the  aorta.  As  seen  in  figure 
2,  both  these  valves  have  three  pouchlike  leaflet  sections  which,  when 
opened,  press  against  the  arterial  wall.  When  blood  attempts  a rover si 
flow  the  pouches  immediately  fill  and  the  ensuing  blood  pressure  for  • s 
losure  ot  the  three  leaflets  and  prevents  backflow. 


From  Structure  Ami  Function 
nt  The  Heart. 

Copyright  1 1964  by  bo'.  Heath 
and  Co. 
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Ibow  1 1 .sr' i • i 1 1 ■ ■ i!  i vmI  v<  ■ , out  standing  examples  of  natur  * 
engim  ring  ipabilities  that  they  are,  may  be  altered  or  destroyed, 

■revei  , by  rngenital  mal formation  or  by  disease . Generally,  it  is  the 
mitral  or  the  aortic  valve  that  is  affected.  There  may  be  insufficiency 
(regurgitation)  due  to  stiffening  and  shortening  of  the  valve  leaflets 
which  prevents  full  closure,  or  there  may  be  stenosis  due  to  narrowing 
if  the  valve  orifice. 

A description  of  the  complicated  fluid  flow  phenomena  occurring 
through  the  aortic  valve,  as  an  example,  impresses  one  with  the  difficult 
problem  proposed  for  the  computer  graphics  investigation.  Valvular  and 
vascular  flows  encounter  eddying,  vortex  flow,  stagnation  regions,  jet 
impingement,  both  laminar  and  turbulent  flow,  and  even  cavitation.  If 
one  models  the  movement  of  an  aortic  valve  even  in  relatively  a simple 
fashion,  it  is  seen  that  at  low  rates  of  pulsatile  flow,  a flow  trail 
forms  from  the  valve  edge  at  either  side  of  the  valve  orifice,  with  flow 
becoming  turbulent  downstream.  With  an  advancing  flow  rate,  a central 
turbulent  stream  is  encountered  through  the  open  valve,  extending  into 
the  aortic  arch  and  descending  aorta.  If  the  flow  narrows  in  the  aorta 
because  of  a stenosis,  then  a jet  is  a resulting  possible  formation  at 
the  constriction.  As  flow  rate  again  increases,  a vortex  forms  where 
the  jet  edges  meet  the  aortic  wall.  Continuing  the  flow  increase  nat- 
urally extends  the  jet  downstream  with  increased  turbulence  aiding  in 
changing  the  vortex  direction  in  the  ascendinq  aorta.  The  center  core  of 
the  jet  strikes  the  opening  of  the  innominate  artery,  allowing  two  flow 
branches  to  form  with  one  compounding  the  vortex  motion  and  the  other 
continuing  along  the  vascular  wall.  A stream  of  vortices  forms  alon j the 
aortic  arch  while  eddying  motions  appear  at  branches  downstream  of  t fi< 
aorta.  No  viscous  clamping  of  turbulence  occurs  in  the  descending  mu 
since  flow  velocity  remains  high. 
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Tn  vivo,  these  eddies  that  stream  from  the  valve  edges  are  func- 
tional in  that  they  prevent  the  valve  cusps  from  completely  retracting, 
thus  allowing  the  valve  to  close  quickly  when  required,  and  to  keep 
coronary  vessels  open.  Eddy  trails  are  simultaneously  generated  in  the 
tri-cusp  semilunar  valve  and  are  deflected  in  the  aortic  sinus  and  arch; 
each  contributing  to  continued  eddy  and  pulsatile  flow.  There  is 
evidence  that  the  core  of  the  jet  flow  is  deflected  with  a periodical 
lateral  motion  at  the  wall  due  to  the  wall's  elasticity.  The  advent 
of  helical  flow  motion  is  seen,  generated  by  centrifugal  force  in  the 
curved  aortic  arch  which  subsequently  allows  secondary  flow,  forcing 
turbulent  action  down  to  the  origin  of  main  vascular  branches.  Pulsatile 
flow  with  its  inclusion  of  rapid  minimum-maximum-minimum  velocities 
incurs  brief  backflows.  Turbulence  and  eddy  formation  from  these  sources 
all  add  to  the  description  of  a flow  pattern. 

The  realization  that  surgery  is  not  the  answer  to  salvaging  many 
irreparably  damaged  heart  valves  led  investigators  to  the  thought  of 
replacement  by  artificial  heart  valves.  With  the  advent  of  the  artificial 
heart  valve  it  is  to  be  noted,  however,  that  the  many  fluid  dynamic  pro- 
blems projected  in  a cursory  way  in  the  above  paragraph  are  still  with  us, 
and  the  resolution  of  such  problems  do  not  enter  into  the  design  of  such 
an  artificial  device.  New  models  are  continually  being  introduced;  even 
cadaver  valves  are  placed  in  metal  rings  for  insertion,  yet  dissatisfaction 
is  still  indicated  with  presently  available  valves.  As  Braunwald  and 
Detmer  ^ pointed  out,  no  currently  available  artificial  heart  valve  ful- 
fills all  the  acceptance  criteria  for  an  ideal  substitute,  as  proposed 
2 

by  Harken  , i .e . , 

1.  It  must  not  propagate  emboli. 

2.  It  must  be  chemically  inert,  and  not  damage  blood  elements. 


3.  It  must  offer  no  resistance  to  physiologic  flows. 

4.  It  must  close  promptly  (less  than  0.05  sec.) 

5.  It  must  remain  closed  during  the  appropriate  phase  of  the 

t 

cardiac  cycle. 

6.  It  must  have  lasting  physical  and  geometric  features. 

7.  It  must  be  inserted  in  a physiologic  site,  generally  the 
normal  anatomic  site. 

B.  It  must  be  capable  of  permanent  fixation. 

9.  It  must  not  annoy  the  patient. 

10.  It  must  be  technically  practical  to  insert. 

In  fact,  true  functional  evaluation  and  comparison  of  present 
artificial  heart  valves  is  rare  among  those  who  manufacture  them, 
beyond  the  relatively  unsophisticated  quality  control  procedures. 

This  is  not  to  say,  however,  that  in  general,  a substitute  valve  does 
not  serve  its  purpose  well.  The  patient  has  been  given  a considerably 
improved  outlook  over  the  obvious  alternative,  with  even  the  presently 


available  prosthetic  heart  valves. 


CHAPTER  Til 


PROBLEM  AREA  I:  ARTIFICIAL  HEART  VALVES . 

The  first  successful  implantation  of  a prosthetic  heart  valve  appears 

3 

to  have  been  performed  by  Hufnagel  who  implanted  a form  of  what  ha. 
become  known  as  the  ball -valve  type.  The  success  was  not  spectacular 
since  improvement  to  upper  extremity  arterial  blood  flow  was  not  accom- 
plished. This  was  due  to  the  requirement  of  cross  clamping  of  certain 

arteries,  for  the  heart  lung  machine  was  not  to  be  invented  until  the 

4 

following  year.  It  remained  for  Starr  , however,  to  achieve  impressive 
results  after  Harken  3 implanted  a prosthetic  valve  with  aid  of  the  heart 
lung  machine  in  1960.  Since  then,  some  60,000  artificial  heart  valves 
have  been  implanted  in  Americans  with  varying  degrees  of  success.  Sub- 
sequently, some  50  to  75  types  of  artificial  heart  valves  have  been 
devised,  with  the  Starr-Edwards  ball  type  valve  being  used  most  extensively. 
Still  present  are  problems  of  thromboembolus , red  cell  destruction, 
material  fatigue  and  chemical  change,  leakage,  and  infection. 

Mention  has  been  made  of  a "ball-type"  prosthetic  device.  Basically, 
prosthetic  heart  valves  may  be  classified  as  those  with  central  flow  or 
lateral  flow.  The  natural  leaflet  valve,  as  shown  in  figure  3c,  allows 
an  unobstructed  central  pathway  for  blood  flow  as  the  leaflets  open. 

On  a cross-sectional  view  an  artificial  leaflet  valve  appears  as  shown 
in  figure  4.  In  attempting  to  replace  the  natural  valve  invest iqators 
found  that  no  man  made  material  can  withstand  the  constant  flexing, 
therefore,  heterograft  (from  other  species)  or  homoqraft  (from  human 
cadavers)  leaflet  valves  are  sometimes  employed.  The  attribute  of  a 
homograft  valve  is  that  there  appears  to  be  a lesser  amount  of  throm- 
boembolism occurrence  than  for  a mechanical  valve,  however,  i statist!  ul 


i 
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. ' i!y  >:  subsequent  lurability  is  still  lacking.  Aortic  insufficiency , 

•i  ;oot  blood  flow  when  replacement  is  made  for  the  aortic  valve  with 
at  . 1 1 1'  aflots,  is  present  due  to  the  difficult  operative  technique 

! p.n  ! • 1,  plus  insufficient  methods  of  preservation  for  the  leaflets. 


1 i.  1,  i , b , Ball,  Di sc (Hammer  smith)  and  Leaflet  Valves. 

There  is  some  evidence  by  Malm  et  al.  J that  leaflets  exposed  to  a high 
• m t qy  • 1<  tron  beam  allows  superior  tensile  strength,  however,  such  a 
>te-ri  1 ization  technique  is  not  widely  available  nor  may  it  be  permanent. 
A problem  relative  to  leaflet  valve  design  is  its  pliability  and  the 
requirement  of  integrated  motion  of  the  leaflets  with  the  bloodstream 
ve lo  l ty  for  re  gating  a shearing  effect  . Consequently,  single  leaflet 
or  multiple  cusp  designs  need  tn  be  precisely  angled  to  the  blood  Mow 


velocity,  and  this  is  presently  not  attainable.  .audits  employing  pals  • 


duplicators  show  one  leaflet  of  a tricuspid  aortic  valve  flexing  with 
the  other  leaflets  allowing  turbulence  and  cavitation  to  be  initiated 
on  the  conforming  side  of  the  valve. 


Fig.  4:  Flexible  Cusp  Type  Artificial  Heart  Valves 

Viewing  the  lateral  flow  classification  of  prosthetic  heart  valves, 
as  mention'd,  the  most  popular  has  been  the  typo  with  a aged  ball,  tor 
the  aortic  position,  as  shown  in  figures  3a  tnd  5,  and  as  discussed  Lat<  r, 
the  caged  disk  or  hinged  form  for  the  mitral  replacement,  as  shown  in 
figures  3b  and  6.  (Although  the  same  typi  of  valve  has  been  greatly  use  1 
for  all  valvular  positions,  in  iiffcrent  sixes,  it:  appears  that  different 
blood  flow  regimes  are  present  in  the  sort  i . mi  mitral  ' nations.  This 
comment  requires  further  investigation).  Interestingly  enough,  as  late 
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a.'  1 Duvoi  sin  md  M Goon  ot  the  Mayo  clinic  mentioned  that  the 

1 1 an  lu  ii  to  .i  ball  type  valve  from  the  difficult  to  manufacture 
1 ■ i fit  • ty»  ••  valve  was  an  easy  decision  to  make  "by  the  real  i ::  it  ion 
that  ?'•  •!  ertain  other  mechanical  functions  found  in  nature,  man  has 
: 1 i 1 ■ : . impn  sm  nts  in  design , su<  !,  is  ' h(  jet  engint  iv<  i the 

• l ipj  in:  w ing  ' tli”  iitd,  and  the  screw  ov<  a th<  fish's  tail".  Rathe i 
v i -s  Lified  statement  perhaps,  since  no  mention  was  made  of  fluid 
lynami  ncepts  required  for  optimized  flow  and  the  still  unformed 
tud  • i t in  ait: nut e lesigned  prostheti  levice . 


5:  Diagram  Showing  Flow 

Direction  Through  A 
Ha l 1 Type  Art  i f iei a l 
Heart  Valve. 


Fig.  <i : Diagram  Showing  Flow 

Direction  Through  A 
Disc-Type  Artificial 
Heart  Valve. 


Fi  i. 
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This  most  widely  used  artificial  valve,  the  Starr-Edwards  typo, 
has  a silicon  rubber  (silastic)  ball  which  is  forced  away  from  the  valve 
ring  position  for  desired  flow  and  which  is  seated  against  the  valve  ring 
to  prevent  reverse  flow  at  the  proper  times.  The  cage  of  this  valve  is 
bulky,  however,  causing  a problem  with  mitral  valve  replacement  for  some 
patients  because  of  their  left  ventricular  chamber  contour  when  in  the 
collapsed  position.  The  cage  may  cause  septal  endocardial  irritation 
and  cardiac  arrhythmias  while  the  ball  could  interfere  in  the  left 
ventricular  outflow.  In  turn,  the  valve,  in  its  entirety,  may  interfere 
with  myocardial  contraction  and  be  related  to  an  ineffective  cardiac 
output.  Maximum  hemodynamic  benefit  is  probably  nullified  and  throm- 
boembolism may  be  realized  due  to  valvular  obstruction.  Post-operative 
autopsies  have  shown  linear  myocardial  rupture  along  the  posterolateral 
strut  of  the  valve  cage  with  endocardial  grooving  present  in  the  vicinity 
of  other  cage  struts.  Also,  in  mitral  regurgitation  patients  the  ven- 
tricular septum  mav  be  hypertrophied.  Since  this  would  cause  shortening 
of  the  distance  between  the  mitral  annulus  and  the  septum  in  the  outflow 
tract  region,  the  conduction  system  may  contact  the  artificial  valve's 
age  which  is  directed  towards  the  septum,  with  trauma  resulting.  Too, 
in  all  valves,  an  interface  traps  red  blood  cells  when  in  the  closing 
phase.  To  prevent  hemolysis,  natural  valve  components  allow  broad 
surface  contact  between  pliable  elastic  tissues.  Hemolysis,  however, 
is  aided  by  a rigid  ball -type  valve  diaphragm  striking  a metal  seat. 

It  is  also  possible  for  the  ball  to  rotate  within  the  valve  cage, 
possibly  damaging  blood  cells,  although  some  investigators  feel  that 
large,  shearing  stress,  present  between  flu.'  ball  and  cage  struts,  allows 


olood  damage  to  occur. 


• i!  variant:*'  in  the  Star r-Kd wards  valve  lias  led  to  the  replacement 
: .....  i ball  in  mi  ode]  by  a hoi]  >w  mi stal  i >ppet , mad 
. i | y railed  Stellite*,  not<  i I i high  » t 1 1 ■ jth  ind  con  i n 
: s*  m . Also,  tin-  -age  struts  and  the  scat  of  the  valve  have  now  beer. 
.*•  r d by  Teflon  cloth  in  an  attempt  to  form  tissue  ingrowth  ind  prevent 
• . . ;■  ition.  rhi  1 ng-t<  rm  effect  oi  tissue  ingrowth  is  unknown 
chr  ml  eml  mi  1 : r si  nt.  it  is  noted  that  some  experiment  - 

,n  is  being  performed  at  i few  installations  with  vitrous  carbon  r 
i ; a less  steel  for  ball  replacement. 

Q 

As  t refine of  thi  ball  valve  design , Mag  >vi  m li  Lgned  a 

ithi  ■■  th  multipli  pin  fixation , as  shown  in  figure  7 . 


!.  7:  Maqovern  Artificial  Heart  Valve 


t,  ||ihn  m ,i|  |..y  oriristing  of  Cobalt  (approx,  hit).  Chromium 
( . , . , '4  d ybd*  ti um  ( ')  -i  ■ ) , r 4 i * k * 1 ( 1 . 75-J  . 7‘>t)  . 
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The  pins  replace  sutures  at  the  valve  ring  position  since  sutures  are 

time-consuming,  tedious,  and  prolong  the  duration  of  bypass  of  blood 

circulation  through  the  heart-lung  machine  during  surgery.  Although 

success  with  this  valve  is  comparable  to  the  Starr-Edwards  type,  it  is 

less  widely  used  due  to  incidences  of  infection,  detachment  of  pins, 

and  bending  of  pins  if  calcified  deposits  are  present  at  the  valve 

seating  position.  Another  refinement  concept  was  sought  by  development 

9 

of  the  Smelof f-Cutter  Full  Flow  Orifice  valve  . In  this  design,  as 

seen  in  figure  8,  a double-caged  valve  was  designed  to  increase  the 

ori fice-to-bal 1 diameter  ratio  so  as  to  allow  more  of  an  orifice  for  flow. 


: Smo  lo  f f -Cut  ter  Full  Flow  Orifice  Artificial  Heart  Vulv< 


Fig.  H 


lb 


Sine  the  diameter  of  the  ball  is  similar  to  the  diameter  of  the  orifice, 
obstruction  to  blood  flow  across  such  a valve  was  assumed  to  be  minimal, 
"hi  c»i':  1 a ranee  of  the  second  caqe  upstream  of  the  orifice  allows  stoppage 
of  motion  of  the  ball  rather  than  having  the  sphere  seat  on  the  orifice 
i in.  Too,  open-ended  cages  are  intended  to  prevent  areas  of  stagnation 
will  n ran  lead  to  embolism  formation.  In  practice  a disadvantage  can 
■ Mir  with  such  a design  that,  is  not  present  in  the  Starr-Edwards  config- 
uration. If  a small  change  occurs  in  the  configuration  of  the  ball,  the 
ball  can  impact  in  the  inflow  orifice.  Clinical  trials  have  not  been  as 
successful  as  in  vitro  experiments. 

With  the  advent  of  the  ball-type  prosthetic  heart  valve  also  came 
the  realization  that  ball  valve  cages  are  bulky,  causing  a problem  in 
some  patients  because  of  the  contour  of  the  left  ventricular  chamber, 
besides  the  aforementioned  hypertrophy  problem.  A low  profile  valve 
is  shown  in  figure  6.  One  variance,  however,  is  an  angled  one  such  as 
the  hammersmith  type  previously  shown  in  figure  3b.  At  the  present 

tine  there  are  perhaps  a score  of  such  disc  valves;  all  with  low  profile 

10 

cages.  Again,  it  was  Hufnagel  that  proposed  the  first  valve  of  this 

nature  in  1965.  This  valve  included  a polypropylene  frame  and  silicone 

coated  disc.  Other  types  are  the  Cross-Jones  valve  ^ with  an  open- 

ended  cage  and  silicone  rubber  lens,  reinforced  with  a titanium  inner 

ring  (figure  9),  the  Kay-Shiley  valve  with  two  parallel  struts  and  a 

plain  silicone  rubber  disc  (figure  10),  the  Me  1 rose-Al varez  valve  ' ^ 

14 

with  a polypropy lerx  ring  and  occluder  (figure  11),  the  Gott  valve 
with  heparin-benzyl -coated  silicon*  hinge  flaps  (figure  12)  and  the 
Davile  valve  ' with  felt  ring  and  polypropylene  occluder;  figure  13. 


Cross  Jones 
Artificial 
Heart  Valve 


Fig.  10:  Kay-Shiley 
Artifi  ial 
Heart  Valve 


Fig.  11:  Me lrose -Alvarez 
Art i ficial 
Heart  Valve 


Fig.  I i:  Davila  Artificial 
Heart  Valve 


<«>»  t Art  i f ici  al 
Heart  Valve 
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Fiq.  16:  Cooley-Bloodwe 1 1-Cutter  Artificial  Heart  Valve 

To  compound  the  problem  of  proper  choice,  there  are  also  design 

variations  such  as  the  hinged  leaflet  shown  in  figure  17.  Two  such 

16 

examples  of  this  style  are  the  Li  1 lehei-Kaster  Pivoting  disc  , 
figure  18,  and  the  Wada-Cutter  valve  ' ^ , figure  19.  Personal  experience 


has  shown  flow  defects  in  the  Wada  valve. 
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Son* c study  has  also  been  accomplished  employing  a true  poppet  type 

valve  as  shown  in  figures  20  and  21.  The  Barnard  UCT  A06  valve  is  of 

1 B 

this  type.  Vorhauer  and  McElhaney  developed  a tear-drop  shaped 
artificial  heart  valve  which  has  some  merit  due  to  its  streamlining 
for  aid  to  blood  flow,  however,  until  the  original  design  is  amplified, 
the  design  performs  a disfavor  to  a collapsed  ventricle  by  the  length 
and  sharpness  of  the  occluder  and  protective  vanes.  In  fairness,  it 
should  be  mentioned  that  the  developers  propose  it  only  as  an  aortic- 
valve  replacement. 


Fig.  20:  Diagram  Showing  Flow  Fig.  21:  Barnard  UCT 
Direction  Through  A Artificial 

Poppet -Type  Artificial  Heart  Valve 

Heart  Valve 


L. 
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Fiy.  22:  Vorhauer  Artificial  Heart  Valve 

A final  example  of  the  many  desi  in  variations  of  the  pros' heti 
heart  valve  is  the  Nakib  toroidal  heart  valve  ' as  shown  in  figure 
23.  this  is  an  interesting  design  since  it  includes  a titanium  ring 
t.i  it  is  pushed  away  from  the  centered  and  stationary  poppet  upon 
ng.  his  allows  cent ra 1 ind  i Iso  peripheral  flow,  yet  ; 

• 1 1 ■ i . : • • • : Lt(  lir<  tion  blo< "l  flow.  Th< ire  is  ilso  a • • lu  tl 

m M-cluder  weight  in  comparison  to  a comparable  ball-type  valve.  To 
: ; ■ • , , imewhat  streamlined  foi  prope j flow  effects . it  has 

in  open-end  low  profile  cage,  however,  lift  1<  follow-up  of  patient 
implantation  is  available  and  no  analytic  studies  have  been  seen  by 
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Fig.  23:  Nakib  Toroidal  Artificial  Heart  Valvo 

Filially,  one  notes  that  the  arguments  put  forth  for  either  the  ball 

or  disc  type  designs  appear  strong  for  either  side.  It  has  been  suggested 

that  lens  (disc)  wear  is  not  as  evenly  distributed  as  ball  wear  under 

20 

comparable  stress  situations  . Since  a ball  tends  to  spin  during  each 
cycle,  the  full  impact  of  the  stress  is  both  reduced  and  distributed 
more  evenly.  In  addition,  a disc  tends  to  move  up  and  down  more  often 
than  it  rotates,  and  also  strikes  the  frame  harder  during  each  cycle. 

At  present,  these  are  theoretical  objections  which  are  not  verified  by 
clinical  tests  us  yet  . Interestingly,  uneven  wear  of  the  disc  has  not. 
boon  shown  cl  ini  cully  when  this  valve  type  is  in  the  mitral  position, 


however,  degeneration  of  the  disc,  when  made  of  a silicon,  is  seen  when 


if  i in  t;u-  lort  ic  position.  On  the  other  hand,  some  studies  have  shown 
a funct  ionu  1 superiority  of  the  disc  valve  over  the  ball  valve,  especially 
in  tin  mitral  section.  The  disc  valve  opens  and  closes  more  quickly, 
then  i s less  regurgitation  during  the  closing  phase  and  the  lens  strikes 
• valvt  si  at  with  Less  force , which  may  negate  a previous  theoret  i i] 
argument.  It  also  has  a shorter  traveling  distance  than  the  ball  valve. 

notes  that  in  the  ball  valve,  chemical  analysis  has  shown  that 
tor  the  implantation  time  of  some  duration,  a few  valves  have  undergone 
bo-  hetniral  change  which  can  then  result  in  mechanical  change.  Such  i 
variant  ball  has  revealed  the  presence  of  free  fatty  acids,  glyceryl 
fatty  acid  esters,  methyl  esters  of  fatty  acids,  cholesterol  fatty  a id 
esters,  and  palmitic  acid.  Also,  rosette  type  crystals,  not  a by-product 
of  vulcanization,  have  been  found.  It  has  also  been  shown  by  the  leading 
manufacturer  of  this  type  of  heart  valve  that  an  affected  ball  shows 
cracks  or  fissures  well  below  the  surface  of  the  ball.  It  is  possible 
then  that  the  biological  change  occurs  with  subsequent  deposit  of  lipid 
material  and  crystals,  and  a large  resulting  deep  fissure  forms;  the 
latter  then  spreads  to  the  surface  of  the  ball.  A requirement  appears 
to  be  then  to  study  the  stresses  as  formed  by  the  action  of  the  arti fi- 
ll 1 heart  valve  components  with  a view  towards  minimizing  such  stresses 
and  the  subsequent  results  when  designing  the  optimum  valve.  Surprinsingl 
perhaps,  and  certainly  worrisome  is  the  fact  that  lethal  degeneration  of 
tht  silastic  ball  in  a Starr-Edwards  prosthetic  valve  type  is  now  coming 
t light  in  some  patients  after  an  implantation  time  of  six  and  one  half 
y-  irs  *.  Figure  24  shows  deformed  and  fractures  poppets  which  either 
remain  immobile  within  the  cage  or  split  into  fragments.  In  passing 
it  is  noted  that  the  senior  writer  of  this  papier  has  incorporated  ex- 
perimental phot ost ross  studies  in  conjunction  with  the  computer  graphics 
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studies  so  as  to  gather  further  of  concern  to  stress  problems,  amongst 
other  requirements. 

In  reiteration,  considerable  confusion  exists  in  the  minds  of 
surgeons  as  to  which  heart  valve  substitute  will  eventually  prove  most 
pre  ferable . 


Fig.  24:  Lethal  Mitral  Valve  Variance 


CHAPTER  IV 
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PROBLEM  AREA  II:  ATHEROSCLEROSIS 

Arteriosclerosis  is  the  term  given  to  a chronic  disease  characterized 

by  abnormal  thickening  and  hardening  of  the  arterial  walls.  Our  interest, 

however,  is  in  atherosclerosis  which  is  a particular  form  of  arteriosclerosis 

formed  by  the  deposition  of  fatty  substances  in  and  fibrosis  of  the  inner 

layers  of  the  arteries.  Since  the  presence  of  atherosclerotic  conditions 

is  one  of  the  prime  causes  of  heart  disorder  and  circulatory  problems, 

then  this  fact,  plus  the  knowledge  that  500,000  people  in  the  United 

States  alone  die  of  heart  attacks  (plus  200,000  of  strokes),  convinces 

one  that  intensified  studies  of  this  disease  are  indeed  important.  Coup- 

22,23 

ling  this  conviction  with  the  thought  that  a number  of  investigators  ' 
suggest  a relationship  between  hemodynamic  factors  and  atherosclerosis 
development  allows  one  to  form  a fluid  dynamic-computer  graphic  inves- 
tigation. 

It  is  startling  to  realize  that  laboratory  studies  with  experimental 

animals  and  postmortem  studies  of  human  infants  have  shown  that  develop- 

24 

ment  of  atherosclerosis  often  appears  shortly  after  birth  . Extensive 
formations  of  atherosclerotic  plaques  were  found  in  American  casualties 
of  the  Korean  war,  as  young  as  a median  age  of  23  years.  At  the  other 
end  of  a history  spectrum,  typical  atherosclerotic  lesions  have  been 
found  in  the  aortas  of  3500  year  old  Egyptian  mummies. 

Campbell  Moses  ’ has  reviewed  certain  historical  investigations 
that  lend  agreement  to  the  thesis  that  mechanical  and  physical  factors 
may  play  roles  as  primary  causal  i ve  I actors  in  the  development  of  athero- 
sclerosis. As  examples,  Carrel  and  Guthrie  in  1906  attempted  to  relate 
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arterial  pressure  to  acceleration  of  atherosclerosis,  while  Klotz 


in 
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1 1 ■■■■•  studied  uncompli  ated  hypertension  and  described  hyperplasia  of  both 

the  intima  and  muscular  elements  of  the  studied  media.  In  1914, 

28 

Anitschkow  discussed  the  accelerating  role  of  hypertension  concerning 

experimental  atherosclerosis,  in  the  presence  of  hypercholesterolemia. 

29 

Allbutt  in  1915  went  further  by  claiming  that  hypertension  is  a prime 
factor  in  atherosclerosis.  Beginning  to  somewhat  employ  mechanics, 

Duguid  1 in  1926  discussed  lesion  localization  at  the  entrance  to  branch 
arteries  when  shearing  stress  occurred  in  the  main  vascular  channel.  More 
recently,  Willis  investigated  localizing  factors  in  atherosclerosis  with 
conclusions  that  hydrostatic  pressure  and/or  lumen  radius  and  curvature  are 
favorable  factors  for  local  atherosclerotic  formations. 


Fig.  25:  Diseased  Atherosclerotic  Human  Artery  (Below)  As  Compared 
To  A Healthy  One  (Above) 


Of  closer  interest  to  our  own  approach  are  the  continuing  studies  of 
32-35 

Texon  . His  hemodynamic  investigations  inevitably  point  toward  the 

laws  of  fluid  mechanics  as  the  prime  "raisons  d'etre"  for  localization  of 
atherosclerotic  lesions  at  specific  areas  of  predilection  in  the  arterial 
tree  as  well  as  for  their  absence  in  identical  adjacent  areas  of  intima. 
Such  areas  of  predilection  are  specifically  declared  to  be  sites  of  dimi- 
nished lateral  or  static  pressure  which  occur  in  arteries  at  sites  of 
tapering,  bifurcation,  branching,  curvature  an!  external  attachment. 
Combinations  and  geometric  variances  of  such  areas  add  to  the  sum  of 
such  patterns.  With  such  a localized  decrease  in  static  pressur  at  a 
site,  occurence  of  a suction  effect  upon  the  wall  is  seen,  with  sub- 
sequent subjection  of  the  arterial  intima  to  the  lifting  effect  of  blood 
flowing  past  upon  the  endothelial  layer  and  subjacent  cells.  Initially, 
the  response  is  a reactive  thickening  by  proliferation  of  endothelial 
cells  and  fibroblasts  from  subjacent  layers.  Successively,  continuing 
blood  flow  enlarges  such  an  initial  plaque,  lateral  pressure  has  continual 
dimunition,  and  cellular  elements  and  lipids  are  added  to  successive 
layers  of  fibroblastic  proliferation.  When  such  pathologic  degeneration 
reaches  a critical  point,  portions  of  the  atherosclerotic  plaques  may 

J 

detach  from  the  arterial  wall,  allowing  embolization  an  1 thrombosis.  With 
such  detachment  the  moving  portion  of  a plaque  may  travel  until  it  reaches 
a section  of  a smaller  artery  and  plugs  it.  If  this  plugging  occurs  in 
the  coronary  artery  it  produces  a heart  attack.  Occurrence  in  the  brain 
causes  a cerebral  stroke  while  gangrene  is  the  result  of  blockage  in  an 
artery  ol  the  lower  extremities. 

In  the  .above  short  expose  there  has  not  been  the  attempt  to  assume 
that  the  bases  of  sites  for  predilection  for  atherosclerotic  plaques  are 
completely  explained  by  fluid  mechanics  concepts.  There  are  many  medical 
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investigators  who  accept  as  conclusive  the  reported  evidence  of  a rela- 
tionship between  atherosclerosis  and  elevated  serum  cholesterol  levels. 
The  possible  consequences  of  diet  and  hereditary  factors  concurrent  with 
elevated  serum  cholesterol  levels  are  duly  noted,  however,  our  interest 
and  capabilities  lie  with  fluid  mechanics  arguments  for  atherosclerotic 
pathogenesis . 


CHAPTER  V 


MATHEMATICAL  APPROACH 
A.  Some  Preliminary  Fluid  Flow  Concepts 

/ 

In  initiating  hemodynamic  studies  of  the  circulation  a useful  pre- 
liminary concept  is  to  assume  a long  straight  tube  with  a steady  flow  of 
fluid  inside.  Maintaining  this  steady  flqw  is  a constant  pressure  which 
is  required  because  of  the  fluid's  inherent  viscosity.  Poiseuille,  a 
physician,  in  1842,  was  the  author  of  an  experimental  law  governing  the  flo 
of  such  an  encased  fluid.  Poiseuille  stated  that  the  head  of  the  required 
pressure  is  directly  proportional  to  pipe  length,  flow  rat-e  and  viscosity, 
and  inversely  proportional  to  the  fourth  power  of  the  radius.  (Girard  in 
1813  had  assumed  a third  power  law  which  was  subsequently  corrected  by 
Hagen  in  1839.)  Upon  the  injection  of  d^e  into  the  liquid  for  a visual 
aid,  it  was  noted  that  along  the  pipe  axis,  liquid  motion  was  quicker  than 
at  the  walls.  The  explanation  was  that  the  fluid  particles  moved  in  a 
series  of  parallel  laminae;  the  lamina  of  fluid  closest  to  the  pipe  wall 
was  stationary  and  each  successive  outward  lamina  slipped  against  the 
viscous  friction  of  its  neighbor,  which  had  progressively  lessened  as  one 
approached  the  pipe  axis.  Such  a group  of  parallel  laminaes  of  fluid  con- 
stitutes laminar  flow  or  streamlines  sin^.-e  particles  of  flow  in  one  such 
ideal  streamline  do  not  become  part  of  another  parallel  streamline  of  flow. 
It  is  of  historical  interest  to  note  that  Poiseuille,  as  a physician,  was 
interested  particularly  in  blood  circulation.  By  choosing  qlass  tubes  of 
capillary  size  for  his  flow  experiments  he  was  fortunate  in  obtaining 
laminar  t low  to  an  accurate  degree  which  would  not  have  no cess  . ly  occur  re 
in  arterial  size  tubes.  He  also  chose  water  as  the  working  tluid  since  he 
could  not  prevent  blood  from  coagulating,  and  again  was  fortunate  for  he 
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thuc  avoided  the  anomalous  viscous  properties  present  in  blood.  It  is  now 
known,  however,  that  i f the  flow  rate-  is  increased,  there  is  a point 
reached  where  the  neat  pressure-flow  relationship  breaks  down  and  Laminar 
flow  ceases.  Fluid  particles  then  follow  random  paths  across  the  tube 
as  well  as  a main  movement  through  the  tube,  with  the  result  that  flow 
becomes  turbulent. 

The  rigid  pipe  concept  is  a preliminary  one  for  arterial  flow  studies. 
The  Poiseuille  law  that  one  usually  encounters  has  been  placed  in  a form 
that  is  amenable  to  Navier's  theoretical  equations  for  viscous  flow,  which 
were  later  amplified  by  Stokes  in  the  1840's.  Poiseuille's  equation  is 
questionable  when  it  is  applied  to  the  vascular  system.  In  the  vicinity  of 
branch  points  the  velocity  profile  may  not  fit  the  one  formed  by  the 
equation.  In  arteries,  the  flow  may  not  be  irrotational , while  in  the 
veins,  blood  may  appear  as  a non-Newtonian  fluid,  wherein  apparent 
viscosity  decreases  with  a linear  velocity  increase  and  a vessel  diameter 
increase.  For  higher  hematocrit  values  there  are  greater  viscosity  changes. 

B.  Determination  of  Flow  About  an  Obstacle. 

Turbulence  is  one  of  the  most  important  aspects  of  the  science  of 
fluid  mechanics  and  fluid  dynamics,  and  one  of  the  least  understood;  a 
phenomenon  that  defies  attempts  at  complete  analysis.  Rarely  is  it  pos- 
sible to  produce  fluid  flow  having  considerable  gross  movement  but  no 
turbulence.  Although  laminar  flow  may  be  desirable  and  achievable,  the 
presence  of  a fine  scale  of  fluctuations,  still  having  turbulent  properties, 
may  be  inherent.  To  compound  the  problem,  the  mathematical  treatment  of 
fluid  motion  has  led  to  a set  of  simultaneous  nonlinear  partial  differential 
equations  of  the  second  order;  the  Navicr-Stokes  equations,  together  with 
the  equation  of  continuity.  These  Navier-Stokes  equations  of  motion  are 
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presently  intractable  in  that  no  mathematical  technique  is  available  for 
treating  the  difficulty  of  nonlinearity.  Known  particular  solutions  are 
those  exact  ones  related  to  certain  limiting  cases  of  very  low  speeds  and 
the  inviscid  fluid  situation  with  irrotational  motion.  However,  modified 
Navier-Stokes  equations  may  be  solved  by  various  methods  employing  boundary 
layer  theory,  in  certain  instances. 

Wind  tunnel  experiments  show  that  a fluid  behaves  in  a series  of 
fairly  regular  changes  when  flowing  around  a fixed  object  in  a stream  of 
motion.  At  low  velocity,  laminar  flow  exists  about  the  object  up  to  a 
point  of  critical  speed  where  the  onslaught  of  turbulence  is  then  noticed. 
Vortices  are  shed  alternately  from  the  rear  top  and  bottom  sections  of 
the  immersed  object,  resulting  in  an  oscillating  wake  familiarly  known  as 
a Von  Karman  vortex  street.  Further  increase  in  flow  velocity  causes  the 
vortices  to  deteriorate  into  random  eddies,  signifying  fully  develop*  i 
turbulence.  Thus,  if  such  a pattern  can  be  developed  by  use  of  a computer 
and  subsequently  plotted  or  displayed  optically,  a calculated  experiment 
may  be  formed. 


1.  Navier-Stokes  Equations  Forms  Required. 


The  forms  of  the  Navier-Stokes  equations  required  for  the  present 
investigation  are  those  that  describe  two  dimensional,  incompressible 
viscous  fluid  flow,  namely; 
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wh»re  is  the  kinematic  viscosity,  and  where  the  x and  y symbols  are  the 
usual  Cartesian  coordinates.  Pressure  is  p and  t is  time.  The  velocity 
components  along  the  x and  y axes  are  u and  v. 

The  equation  of  continuity  is 


t>p  ^ op U + dp V _ 
3t  lx  + 3y 


which  is  equivalent  to 


dp  3u  3v 

-r—  + p t — + r — = 0 

dt  >x  <y 
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(4) 


whe  re 


da 

dt 


3 9 3 

9t  + u 3x  + v 3y 
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These  equations  express  the  conservation  of  mass.  Assuming  the  fluid 

dp 

to  be  incompressible  as  noted,  then  — = 0,  and  the  equation  of  con- 


tinuity reduces  to 


3u  3v 

3^+  37=  ° 
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The  stream  function  ip  is  defined  in  terms  of  the  veloc-Lty  components 


so  that  it  satisfies  the  equation  of  continuity.  Thus 

3 ip 

u = — 

3y 


(7) 


and 


3 p 

v = - t n? 


(8) 


The  vorticity  w is  defined  as 

u,  _ _ _1  I 3u^  3jv  1 (9) 

2 3y  " 3 x I 

If  one  differentiates  equation  (1)  with  respect  to  y and  equation  (2) 

with  respect  to  x,  and  notes  that  p is  constant,  then  substracts,  the 

vorticity  transport  equation  is  formed; 

dm  _ 3m  3m  3 ip  3m  '.lip  _ 2 (10) 

dt  3t  + 3x  3y  3y  tx  V ,1J 
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The  stream  function  is  useful  in  determining  flow  streamlines  and  the  vor- 
ticity  function  is  a function  of  fluid  velocity  in  the  vertical  and  hori- 
zontal axes  directions.  As  seen  by  equations  (7)  and  (8),  the  vorticity 
function  is  usually  transformed  into  stream  function  values  to  ease  the 
solving  of  the  fluid  motion  equations.  With  the  aid  of  the  stream  and 
vorticity  functions  the  newly  formed  vorticity  transport  equation  allows 
one  to  have  a fluid  diffusion  equation  in  terms  of  vorticity. 

An  equation  for  p may  be  obtained  by  differentiating  equation  (1)  with 


respect  to  x and  equation  (2)  with  respect  to  y.  The  equation  formed  is 

(11) 


V2p  = 2p 


r 2 2 2 

1 3 \p  3 <P  3 <p 


3x"  By""  3x3y 


On  a fixed  solid  boundary 
u = v = 0, 

however,  in  terms  of  the  stream  function 

d ip  cljp  _ Q 
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where  s and  n refer  to  the  tangential  and  normal  directions  respectively. 
Employing  equations  (1) , (2) , (6)  and  (9) , one  sees  that  on  the  boundary 
of  a solid 

3w  3p  _ 3u  (14) , (15) 
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If  it  is  now  assumed  that  U is  a constant  velocity  and  L is  a char- 
acteristic length  of  the  system,  then  dimensionless  variables,  denoted  with 

a subscript  d,  are  introduced  as  follows: 
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nematic  viscosity 
be  written  as 
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as  before. 


The  vorticity  transport 


(17) 


and  the  vorticity  equation  becomes 

V^tp  = -2(u  (18) 

d d 

The  subletters  attached  to  the  dimensionless  variables  will  be  dropped 
for  convenience  and  all  equations  will  then  be  assumed  to  be  expressed 
in  terms  of  dimensionless  variables. 

The  boundaries  involve  the  specification  of  the  velocity  of  the 
fluid  at  all  boundaries.  To  see  how  the  boundary  conditions  are  spe- 
cified, the  simple  example  of  fluid  flowing  over  the  surface  of  a flat 
plate  is  noticed  at  this  point,  however,  boundary  conditions  will  be 
more  fully  discussed  later.  Thus, 

u (x,0,t)  = 0 
P 

v (x,0,t)  = V g(t)  f(x)  (19) 

P P 

(where  p stands  for  physical  variables)  are  imposed  on  a boundary  at 
y = 0,  then  the  dimensionless  problem  should  have  the  boundary  conditions 

|^(x,0,t)  = 0 (20) 

dy 

X 

i|)(x,0,t)  = const  + — - g(t)  1 f(a)da  (21) 

v o 

where  t is  a dummy  variable.  The  boundary  conditions  contain  the  dimen- 
sionless term  UL/v.  This  is  usually  seen  as  which  is  the  Reynolds 
number.  The  Reynolds  number  enters  the  problem  through  the  boundary 
conditions.  The  dimensionless  function  g(t)  allows  a velocity  change 
with  time  on  the  boundary.  It  is  usually  seen  in  the  form 
g(t)  = (l-r."at) 

The  dimensionless  function  f(i)  allows  spacial  variation  of  the  velocity 


on  the  boundary. 
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The  initial  conditions  are 

u(x,y,0)  = 0 

v(x,y,0)  = 0.  (22) 

Thus,  in  equations  (17)  and  (18)  one  has  two  simultaneous  second  order 
partial  differential  equations  in  two  unknowns  ip  and  w with  boundary  and 
initial  conditions  of  the  forms  (20,  (21) , and  (22).  Equation  (17)  is 

a nonlinear  equation  which  allows  analytical  solutions  for  only  very  simple 
flow  problems. 

2.  Computational  Mesh 

To  perform  the  desired  flow  calculations  an  obstacle  was  placed  in 
the  center  of  the  flow  stream.  The  area  about  the  obstacle  was  divided 
into  square  cells,  which,  in  their  entirety,  formed  a basic  computing 
mesh.  An  optimum  arrangement  of  cells  was  formed.  The  fixed  object, 
about  which  flow  was  to  be  studied,  was  placed  with  its  boundaries  along 
mesh  lines.  Enough  cell  spacing  was  allowed  upstream  of  the  obstacle  for 
initiation  of  flow  perturbations  and  initial  conditions  as  required.  (See 
figure  26) 

Figure  27  shows  an  enlarged  portion  of  the  computational  mesh. 

All  partial  differential  equations  involved  had  their  derivatives 
replaced  by  finite  difference  equivalents  !\  which  allowed  a series  of 
finite  difference  equations  instead  of  the  difficult  differential  equations. 
These  finite  difference  approximations  are  algebraic  in  form  and  relate 
the  value  of  the  dependent  variables  at  a point  in  the  mesh  to  the  values 
at  a number  of  symmetrically  arranged  neighboring  mesh  points,  at  one 
instant  of  time  and  then  succeeding  time  intervals.  Although  the  method 
sounds  time  consuming,  and  can  be,  if  performed  by  hand  for  approximately 
the  1500  cells  available  in  the  present  mesh  network,  the  method  becomes 
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q'.t't-e  r i>  ini'  and  immensely  accelerated  when  the  method  of  squares  or 

■JQ 

:•  i ixat  i n m thod  is  applifad  with  their  ramifications,  and  in  turn, 
pj  jramm  I foi  the  ligital  >mput  er.  Cal  u 1 at  i or  t at-  t hen  mat  ter:; 

! st  >;u)s.  Values  of  the  velocity  components  were  calculated  as  cen- 
teri  i lifferen  • 1 th<  stream  function  at  each  cell  corner. 

Figure  28  shows  tin  irregular  computing  mesh  employed  so  as  to  allow 
jeneralizt  1 lii  i inal  juantit  < s * ■ f ■ ng  the  difference  equations . 

As  will  be  seen,  this  generalized  concept  was  useful  for  developing  the 

y requii  l to  analyze  obstacles  that  consisted  of  curved  boundaries. 


Fig.  26:  An  Idealized  Obstacle  and  its  Computational  Mesh 
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27;  Portion  of  Computational  Mesh 
o = stream  function. 
u,v  = velocity  components. 
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2H:  The  Irregular  Computational  Mesh 
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3.  Finite  Difference  Equations  Hequired. 


It  was  first  necessary  to  form  finite  difference  approximations  to 
the  vorticity  transport  equation  and  the  equation  linking  vorticity  and 
the  stream  function.  Expanding  by  means  of  a Taylor's  series  formation 
and  referring  to  figure  28  . required  evaluation  about  points  A and  C 
are  seen  to  be 


= uQ  + a|hu_)  + l/2a^h  u^  + 0(b) 

X XX 
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u = uQ  - a3h  u0  t l/2ayUo  + 

X XX 
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where  u is  a dummy  variable  for  ip  or  w. 


Eliminating  u0  one  obtains 
xx 
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by  ignoring  the  higher  orders  of  h. 


Elimination  of  allows 
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In  a similar  manner, 
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Eliminating  u0  one  obtains 

yy 

U0  = + <a2  " a4>  UQ  ‘ a2UD)/ha2a4U2  + V (29) 
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and  similarly 
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Kqn.it ions  (21)  throuqh  ( 30)  .ire  then  used  to  obtain  the  t i n i t • • li  ft.-r.-nf  •• 
approxiraal  ions  to  the  nomlimensional  forms  of  equations  (17)  an  1 ,11'  . 
Placing  equations  (25),  (26),  (29),  and  (30)  in  slightly  lift'  e nt  I : : > 

( 31) 


' 2 2 2 2 
~ = a3UA  + (ara3)U0-aiUC 

3x  q 


ha1 a^ (a^+a^) 


£ 

3x 


Li;  = a3V(a3+ai)  uo+aiuc 


T aia3(ai+a3) 


(32) 


3u 

3y 


aK+(Va4)uo"a2UD 


O 


ha2a4U2  + a4) 


( 33) 


£ 

3y^ 


= a4UB~(a2+a4)Va2UD 


O ^a2a4(Va4) 


(34) 


Recalling  that  the  vorticity  and  stream  function  relationship  is 
Vi 1 = -2w  = 3 ^ if  + 3^iji 

2 J 

3x  3y 

and  the  vorticity  transport  equation  is 

3(p  _ , ,2  + jhji  jho  3jl  3co 

St”'10  "3x  3y"  ~ 3y"  3)t 

or 

2 2 

3oi  _ 3 (o  3 <o  dip  3(o  3i|i  3co 

3t  ” , 2 + 7T  + 7 “ 7~ 

i x 3y  3 x 3 y 3y  3 x 

and  that  the  nodes  A,li,C,D,  and  0 of  each  mesh  cell  have  the  foil  wing 

positions  in  the  i.j  notation 

A=i  + l,j  B - i,j  + l C = i -1  , j D = i , j -1  0=i,| 

then  by  substitution,  (with  the  dummy  variable  u becoming  y and 
in  the  proper  places,  the  vorticity  transport  equation  in  finite  diif  tm  • 


form  is; 


4 2 


Su> 


r am  “ (a  +a  )w.  + a iii, 

? 1+1,3  3 1 1,3  1 i — 1 , 3 


h_  aia3(a1+a3) 


a.u  - (a_+a.)u>.  . + a_oj.  . ■ 

4 1 , 3+l  2 4 1,3  21,3-1 


WW 


2 2 2 2 

a v . , + (a,  -a,)  ifi . . -a,  1J1.  , 

3 1+1,3  l 3 1,3  1 1-1,3 

ha1a3(a)+a?) 


2 . 2 2,  2 

r a4ll,i  , j+l+  (d2_a4)  “i  , j'a2Wi  , j-1 

SVVV 


2 , 2 2,  2 

a3“i+lfj+(ara3,“i,j-ai"i-lfj 

ha2a4(d2+a4)  , 


a4^i,3  + l+(a2-a4)li;i,j-a2^i,3-l  i 


ha,  a (a  +a  ) 
13  13 


Regrouping : 


3u> 

at 


h a^aj+a  ) 


Gl)  — 

1+1,3  . 2 


CO . . + 


h“a3a3  h2a3(a3+a") 


U) . 

3.  i-l, 


! 2 


GJ.  . 


h a2(a2+a  ) 


4i  "i-i+l  u2,  , U’i,  j + . 2 


h a2a4 


h a4(a2+a4) 


00 . . 

1.] 


a3h 


(a3-a3) h 


h a (a  +a  ) 

V 


ip.  , . + 

1 + 1,;)  . 2 


h a3a3 


aih 


1 j “ . 2 


h a3 (a^+a3) 


(35) 


1 i 


V 


h VW 


‘i, j+1 


v (a2~Vh 

^2*4 


U)  . 

1 , 1 


- a2h 


i 2 , % i , j+1 

h a^a^) 


+ a 


u) . 


hVw  i*i’j 


.2 

H aia3 


+ r 3 


aih 

— 

? , v i-1 , j f 

h a3(a1+a3) 


a4h 


h a2U2  + a4) 


, j+1 


+ U2~a4)h 

h2a  a 
2 4 


a0h 


0.  ■ 


h a (a  +a  ) 
4 2 4 


(36) 


2 

In  turn  the  stream  function  V ip  = -2ai  has  the  finite  difference 


a ijj.  . - (a  +a  ) ijj.  . + a,  il>. 

3 i+l,3  3 1 i,j  1 1-1,3 


T aia3(al+a3> 


Vi, j + 1 - (VV  Vj  + Vi,j-1 

h2 

T a2a4(a2+a4) 


-2w . 


1,3 


( 17) 


Regrouping : 
1 


VW  1+1 '3 


1 , + 1 


La3  i i j a3 (ai+a3)  i_1'j 


+ 1 


*4  4 . , ' ~ ^ 4 + 1 


WV  i'3  + 1 a2a4  i,:i  VW  i'3_1 


.2 

-h  uj  . 
i , 3 


(38) 


form 
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■ th<  ora ■:  itei  program  cal  :ulated  the  type  ol  each  grid  point, 

• pi  )!  un  then  calculated  the  mesh  distances  a ,un,a  } and  for  each 
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node.  As  Booy  pointed  out  , the  analysis  of  a difference  equation  on 
points  next  to  a boundary  involves  a finite  difference  stencil  with  arms 
Lc  , l = 1,2,3,4)  of  various  lengths  such  as  occurred  in  our  problem. 

It  was  his  inherent  stipulation  that  no  mesh  distance  a.  be  less  than 
:).25h.  It  was  therefore  decided  that  for  our  program,  if  a mesh  node 
was  within  0.251)  of  the  boundary,  then  that  node  was  redesignated  as  a 
boundary  point  as  shown  in  figure  20. 


Fig.  2l>:  Finite  Difference  Stencil  Near  A Curved  Boundary. 
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4.  Numerical  Procedure 

The  details  of  the  numerical  procedure  can  be  confu.  in  : : n i t 1- 

lowed  closely,  therefore  a generalized  view  will  first  }ivi  , 
by  a detailed  discussion  which  then  culminates  in  the  pi  . at  i t j n : t 

first  set  of  pictorial  results. 

At  the  start  of  a computational  experiment  values  ;i  the  iium  ntal 
and  vertical  components  of  the  average  fluid  velocity,  1 > th.  r with 
vorticity  and  stream  function  values  at  each  cell  nodal  ; int  it  c ■ 
and  at  initial  time  must  be  specified.  Then,  cu,  the  vorti  ity  fui.  -t:  ■ 
and  . , the  stream  function  are  calculated  at  each  corner  of  each  < 
with  the  knowledge  of  initial  conditions  of  the  flow.  The  vorti  ity 
transport  equation,  in  its  equivalent  algebraic  form,  in  th  n > 1 v 1 
so  as  to  replace  initial  vorticity  values  at  all  lattice  point  with 
Ivanced  values  in  successive  times.  Such  are  used  as  source  t<  m ; 
the  Poisson  equation  solved  for  the  stream  function.  The  oqunt i n f i 
the  stream  function  on  the  interior  of  the  flow  field  is  then  solved  i: 
a boundary-value  problem,  with  the  boundary  stream  functi  >n  ; i f i • 1 by 
programmed  instructions  to  the  computer.  Vorticity  values  ir.  ilculat  • 1 
in  an  iterative  fashion  on  the  boundaries  of  the  obstacle  forme  1 within 
the  flow  field  and  the  process  is  repeated  until  the  maximum  hanges  in 
vorticity  and  stream  function  at  any  nodal  point  in  the  field  reach  a 
prescribed  tolerance  of  error  within  the  computational  pro  • . . Aft-  r 

this  tolerance  has  been  reached  the  computer-programmed  pn  • i . - 1 

again  with  new  vorticity  and  stream  function  values  on  the  Ixaunlary.  n 
is  thus  allowed  to  visualize  the  flow  pattern  and  vorticity  -ml  i il«  ' 
the  obstacle  for  changing  increments  of  time.  The  eomput-  r mui.t  , < : 
necessity,  be  supplied  with  operating  instruction-,  for  d ul.it  : nui 
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f ; w >n  figuration  at  any  sequential  instant,  using  the  values  pre- 
vious ly  mentioned  as  data  input  at  each  point.  •Jumerical  stability  must 
t it  ! and  the  requirement  of  truncati  -n  turbulence  not  to  bo  allowed, 
n;  u • The  operating  instructions  are  necessarily  tedious  and  time 
step-by-stej  n tru  ti  ins  in  r<  jui  red , yet  th<  pi  >gram 
is  written  with  sufficient  generality  to  solve  a wide  variety  of  similar 
■ r<  * L-  ns.  For  example,  for  the  heart  valve  study,  an  object  in  the 
•treat  flow  may  be  studied  with  an  approximately  equivalent  shape  to  a 
ie lit  valve  section  of  interest,  but  with  the  allowance  of  a family  of 
like  irdiac-valve  shapes  to  be  written  into  the  computer  program.  Once 
the  verall  format  is  seen,  operating  instructions  in  terms  of  specific 
requirements  can  be  initiated  for  a particularly  shaped  valve  segment, 
but  with  the  allowance  for  changing  requirements  and/or  shaped  segments. 
The  program  will  allow  for  additionai  and  immediate  required  evaluations 
,*  j re scribed  points  in  space  or  time.  Ultimately,  the  results  are  pre— 
nti  i in  any  unount  f tetail  desired  by  means  of  an  optical  lisj  lay, 
iata  print-outs,  or  single  time-frame  plots. 

Detail-wise,  the  numerical  procedure  is  a successor  to  the  technique 
formed  by  Fromm  and  Harlow  11 . Their  procedure  was  basically  as 
fol lows  : 

i)  At  the  beginning  of  a cycle  all  required  variable  values 
were  placed  into  computer  memory  slots  as  initial  lata 
or  results  from  a previous  cycle, 
it)  ince  vorticity  and  stream  function  values  were  assumed 
known  at  each  mesh  point  at  a particular  time,  a Dufort- 
Franke 1 ^ difference  equation  approach  was  utilized  to 
solve  the  vorticity  transport  equation.  To  handle  non- 
linear advection  terms,  a special  difference  form  was 
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• applied  which  conserved  vorticity  and  kinetic  energy,  thus 

preventing  nonlinear  instability. 

$ 

iii)  The  vorticity  function,  in  terms  of  the  stream  funct  i was 

37 

solved  by  relaxation,  using  the  Liebmann  procedur  t 
obtain  new  stream  function  values, 

iv)  The  stream  function  values  were  then  used  to  obtain  new 
velocity  components. 

v)  New  velocity  and  vorticity  values  were  substituted  into 
the  vorticity -transport  equation  and  the  process  was  re- 
peated until  a steady  state  was  obtained. 

In  their  original  studies,  simulated  motion  of  air  past  an  oblong 
obstacle  that  was  involved  in  bomb  blast  studies  produced  a streamline 
plot.  A quite  different  streamline  plot  is  formed  if  the  air  mass,  rather 
than  the  obstacle,  is  stationary  with  respect  to  the  observer.  Thus, 
streamline  and  also  isotherm  plots  are  possible. 

Although  results  of  the  iteration  procedure  were  impressive, 
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Pearson  found  the  wall  boundary  condition  led  to  instability  within 
the  flow  field.  It  also  appears  that  the  Dufort-Frankel  process  for 
solving  the  vorticity  transport  equation  has  no  true  advantage  over 
equivalent  processes.  Computation  time  is  not  minimum  due  to  non-usage 
of  the  over-relaxation  technique.  Also,  the  use  of  periodic  boundary 
conditions  is  artificial. 

44 

Pearson,  and  also  Esch  , of  the  same  laboratory,  formed  improve- 
ments in  the  procedure  by  solving  the  vorticif y-transport  equation  using 

45 

the  Peaceman-Hach ford  method  . Ove r-re laxat ion  was  also  used  to  spec  1 
computations.  Pearson's  and  Ksch's  analyses  of  the  finite-difference 
equation,  excluding  the  nonlinear  adveetion  terms,  have  apparently  shown 
that  the  time  increment  requirement  needs  only  be  not  so  large  a ; to 
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obscure  i he  transients  of  the  solution.  A stability  analysis,  performed 

so  as  to  include  nonlinear  terms,  also  appears  to  bear  out  this  result. 

46 

Weber  of  the  University  of  Utah  agreed  with  Pearson’s  and  L'sch’s  conclu- 
sions concerning  stability  only  so  far  as  the  Reynolds  number  for  flow  was 
small.  Thought  was  given  to  decreasing  the  time  steps  to  prevent  divergence 
at  higher  Reynolds  numbers  (a  similarity  parameter  useful  in  fluid  flow 
studies) . 

Pearson  and  Weber,  to  some  extent,  each  attempted  to  modify  the  iteration 
technique  as  well  as  to  study  the  truncation  error;  the  truncation  error 
being  the  error  inherent  in  the  use  of  integration  formulae,  that  is,  the 
difference  between  tiie  true  and  approximate  solutions  of  a differential 
equation  due  to  one  step  in  the  numerical  integration  process . These 
errors  are  bound  to  accumulate  in  step-by-step  procedures  since  the  starting 
point  at  each  is  approximate  and  each  succeeding  step  compounds  the  error. 
Extensive  calculations  of  concern  to  the  problem  of  truncation  and  of 
concern  to  the  choices  of  smoothing  parameters  and  convergence  criteria 
are  being  studied  in  our  own  investigations.  In  comparing  Weber's  com- 
puter program  modifications  with  Pearson's,  it  appears  that  the  former's 
method  is  more  accurate.  Vorticity  and  stream  function  iteration  type 
errors  appear  to  be  smaller  by  several  factors.  There  may  be  some 
question,  however,  as  to  equal  accuracy  where  higher  Reynolds  numbers 
are  concerned. 

Based  upon  the  above  arguments,  the  approach  of  Pearson,  modified 
by  Esch  and  extended  by  Weber  was  attempted.  In  turn,  programming  sub- 
routines were  added  to  allow  results  for  our  particular  problem,  as  well 
as  to  allow  Gerber  plot  output  and  graphics  console  capability.  The 
solution  space  was  now  divided  into  i computational  mesh  as  shown  in 


f 1 gure  10  . 
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Fiq.  30:  Computational  Mesh  Form  for  the  Pearson-Esch  Approach. 

These  mesh  forms  included  unequal  stencil  arms  as  have  been  seen 
figures  25  and  29  so  as  to  compensate  for  eventual  emplov-.ient  of 
boundaries.  Group  A nodes  consisted  of  boundary  points,  group  B 
tlie  adjacent  points  to  the  boundary,  and  group  C consisted  of  all 
interior  points.  The  mesh  cells  were  declared  to  be 

(x.,y.,t  ) = (ih,jh,k  At)  (39) 

where  \ = 1,2, M 

j = 1,2 M 

k = 1,2 

h Ax  = Ay  = 1 / ( N - 1 ) (40) 


in 

curved 

were 

other 
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discrete  stream  and  vorticity  functions  p.  . and  w.  . were  first 

if j if j 

approximated  by  $(x.,y.,t  ) and  u>(x.,y.,t  ).  It  was  assumed  that  all 

1 J K 1 1 K 

quantities  were  known  at  time  k and  could  be  calculated  for  a time  of 

k+1.  The  procedure  then  became  the  following: 

i)  Superimpose  the  description  of  the  physical  boundaries  on 

tiie  finite-difference  grid  to  determine  the  mesh  lengths 

a^,  i - 1,2, 3, 4,  and  type  the  mesh  nodes  (interior,  adjacent 

k+1 

interior,  and  boundary).  The  values  of  p.  . on  A,  which 

k + 1 

will  be  denoted  as  p , are  taken  to  be  equal  to  the 

specified  boundary  values  of  iK  The  same  is  true  for 
k+1  , .... 

u)^  when  vorticity  is  given  on  the  boundary.  On  the 

k+1 

other  hand,  if  3ij»/3n  is  given  w is  an  unknown  of  the 

A 

problem. 

k+1  k+1 

ii)  Estimate  p and  w by  quadratic  extr  polation  from 

D • v_  A+B+l 


past  values.  The  formulas  are 

k+1  k k-1  k-2 

; . . = 3*.  . - 3<L.  . + ii.  . 

i,3  i»3  i,D  i,l 


(41) 


k+1  k k-1  j k-2 

u) . . = 3ui . - 3<o . . + u . 

i,3  i,3  1,3  i,3 


(42) 


It  is  seen  that  this  involves  storage  of  three  consecutive 

time  levels  of  vorticity  and  stream  function  values. 

k+1 

iii)  Improved  values  of  are  formed  by  employing  the 

formula 


k+1 
ui . 

1,3 


k+ 1 

= smoothing  factor*  u>.  . +( 1-smoothing  factor) 

1,3 


unproved 


•I'1; 


(4  i ) 


2h 


. !k  + 1 

. ,k+l  ,k  + l , k + 1 

-1/2*3,j  + '"l,3  + l+  ♦l.j-l 

_ n ^+1 

- 3h  / 3 p ) 

k + 1 

(44) 

' 9x  1 

| 

1.3 
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The  smoothing  factor  is  necessary  to  prevent  overcorre  tion 
and  subsequent  solution  divergence.  This  formula  is  quickly 
derived  by  taking  an  x=0  and  i=l  boundary  as  an  oxampl'  md 
writing 


, 3 il;  v 2 3 ^ h3  3 3ij/  4 

. = ij<  + h — + h — - + — — - + 0(h  ) 

2,3  3x  . 2 3 

Dx  6 3x 


dtp  4hJ  d2>l>  3h3  3“>  . 4 

+ 2h-2t  + — - — r + ; + 0(h  ) 

3,3  3x  2 _jx2  6 3x3 


(45) 

(46) 


where  all  quantities  on  the  right  are  evaluated  at  (1,  j)  . 

3 3 2 2 

Eliminating  3 ij>/3x  and  then  solving  for  3 y/bx  one 

obtains 


h2  —4  = 4ip  . — 1 / 2 '/.■  . -3.5(1;.  . -3h  — 

3x2  2'3  3'^  1(1  1,3 

1.3 


4 

+0(h  ) 


(47) 


i v) 


by  employing  this  result,  together  with  the  usual  for- 

2 2 2 2 
mula  for  h • 3 ip/dx  in  the  Poisson's  equation  V = -2., 

one  obtains  equation  (43). 

kH 

The  vorticity  transport  equation  is  then  solved  for  • 

D +C 

by  using  the  Poaceman-Rach ford  method.  Recalling  that  the 
working  equations  for  the  stream  function  and  vorti  ity  in 
finite  difference  form  were  shown  in  equations  (36)  md 
(38),  it  is  then  required  in  this  approach  that  tin-  time 
factor  ( 3/3t)  must  be  accounted  for,  at  this  point.  This 
is  accomplished  by  sweeping  across  the  rows  of  nodal  point: 

in  the  mesh  (the  array  is  a matrix  equivalent  as  stored  1 1 ; 

9 

the  computer)  to  solve  for  to.  . at  new  time  values.  For 

1.3 

the  vorticity  transport  equation; 
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Again,  the  procedure  is  followed  for  the  stream  function 

relationship; 

2 I 

V i|>  . = -2ui  , (50) 

,lt  +At  t +At 

I o o 


The  computer  stored  matrix  has  values  of  the  quantity 


uj.  as  defined  by 
if] 


ID  . , r . 5 


k k + 1 i 

US.  . + 0). 

1 3 1.3 


(51) 


where  also 

* 


1.3 


.5 


,k+l  k 1 
ip.  . + ip.  . 

i.3  i . 3 1 


(52) 


One  should  note  that  if  the  stencil  arms  are  equal  for 
each  cell,  as  they  were  for  the  first  set  of  plots  formed 
(a^  i = 1,2, 3, 4 = unity),  then  with  the  application  of 
the  second  central  difference.  62  . 
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the  required  finite  difference  formulae  become 


* k <2  * ,2  k . . 

u) . -a)  o u) . 6 at . . k k k 

1,1  i,i  x i,i  y i,i  1 : w.  , - . . 

— U.  = ~^L  i -y -'-2-  , - ■ i,  j+]  l,  3-l  Vl, 


2 2 2 
. 5At  h h 4h 


I \ 


"'i-l.j  ^i+l.j 


(54) 


k + 1 

0. 

1.3 


* e 2 * r2  k+1 

* o a) . 6 ui . . 

-“i  4 = x *.J  + + I 

1 ' J o 1 2 

h h“  4h 


k + 1 

w . . , 

i,D  + l 


k*  L 


.k+1 
1 1. 

i + 1,  j 


,k  + l 


“ U)  . — (Jj  . j.  - 

1+1,3  1-1,3;  1,3+1  1,3-1 


(55) 


where  i|,.  . uses  the  most  recent  estimate  of  , from 

1,3  i,3 

step  v of  the  procedure.  These  equations  yield  a tri-diagonal 

coefficient  matrix  which  may  be  solved  according  t>  an  al- 
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gonthm  by  Todd 

v)  From  the  nondimensional  form  of  the  vorticity  function, 

k + 1 

improved  values  for  i|;  c are  calculated.  To  accomplish 

this  step  it  is  necessary  to  employ  successive  over- 

, ■ 48  J „ , k+1 

relaxation  and  the  use  of  as  boundary  values. 

It  should  be  noticed  that  the  source  term  in  the  Poisson 

equation  changes,  using  the  results  of  the  previous  step. 

vi)  The  process  is  repeated  by  returning  to  step  iii  and  the 

iterations  continue  until  the  stream  function  and  vorticity 

values  read  a prescribed  maximum  chanqe  toleruno  . The 

tolerance,  as  set  by  Weber,  are  functions  of  the  Reynolds 

number.  Note  that  this  is  because  lie  used  an  absolute  error 

rent  r 'her  t turn  the  more  useful  relative  error  t<  t.  Value 
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of  the  stream  function  and  vorticity  in  the  dimensionless 
forms  in  the  Weber  approach  increase  with  increasing  Rey- 
nolds number.  The  error  may  be  expected  to  increase  also. 

As  to  stability  criteria  it  appeared  to  Pearson  and  Esch  that  no 
riteria  problem  was  presented  as  long  as  the  flow  was  restricted  to 
Reynolds  numbers  less  than  a value  of  100.  Weber's  view  of  stability 
requirements  appear  to  be  correct,  however,  for  higher  Reynolds  numbers. 

For  Reynolds  numbers  in  the  range  of  3000  the  time  step  of  .005  allowed 
convergence  if  the  mesh  size  was  .1.  Recalling  that  the  dimensionless 

function  g(t) , used  to  form  the  velocity  change  with  time  on  the  boun- 
^ ^ 

dary,  was  1-e  , alternately  one  could  restrict 

-20t 

(1-e  ) • (R  - R ) < 30  (56) 

e e 

new  old 

for  non-divergence.  In  turn,  convergence  tolerances  and  truncation 

* 

error  computation  followed  the  approach  of  Weber  . 

5.  Boundary  Conditions. 

Since  the  stream  function  is  specified  as  a boundary  condition, 

: and  either  its  normal  derivative  if  or  the  vorticity  x)  must  also  be 

n 

specified.  The  normal  derivative  was  specified  with  the  stream  function 
at  both  walls  and  at  the  inflow  and  outflow  region  sections  since  phy- 
sically, this  was  equivalent  to  specifying  the  u and  v velocity  com- 
ponents on  the  boundaries. 

* The  investigation  formed  by  Weber  included  the  aid  of  Mr.  Charles  Brauer 
in  forming  an  amount  of  program  subroutines  equivalent  to  certain  numerical 
procedures.  In  turn,  Mr.  Brauer  formed  like  subroutines  for  this  investi- 
gation, thus  allowing  the  Weber  and  Brauer  views  to  be  similar  Ln  both 


projects . 


With 
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the  assumption  of  channel  flow  and  parallel  soli  ! wall:  for  the 

upper  and  lower  boundaries  of  the  solution  space, 

* ( x , 0)  = 0 

* (x,l)  = Ref(t) 

ij>  ( x , 1 ) = tjj  (x,0)  = u = 0 

y 'y  wall 

ip  (x,0)  = rf,  (x,l)  = -v  , = 0 (57) 

x ’x  wall 

i.e.,  the  upper  and  lower  walls  take  on  the  value  of  a stream  li.r  tnd 
the  velocity  at  the  walls  is  zero.  Along  the  left  or  inflow  fcx  iiviary  > 
choice  is  available  for  the  desired  entrance  velocity  profile.  It  is  to 
be  noted  that  rather  than  form  the  normal  derivative  at  the  inflow  and 
outflow  boundaries  as  first  proposed  by  Pearson,  Esch  utilized  the  vor- 
ticity  with  the  stream  function.  In  turn,  by  specifying  the  vorticities 
to  be  equal  to  zero  at  the  entrance  channel,  the  generation  of  unrealisti 
vorticites  at  that  position,  as  is  sometimes  seen  with  the  Pearson  ap- 
proach, was  avoided.  However,  one  could  question  this  point.  Hence, 
the  entrance  profile 

ip  (0,y)  = yR  0 < y < 1 
e — — 

and  thus  i p (0,y)  = u(0,y)  = R_ 

y ^ 

ip  (0,y)  = -v(0,y)  = 0 (58) 

x 

i.e.,  at  first  a rectangular  velocity  distribution  was  specified  at  the 
entrance  to  the  channel.  However,  the  velocity  profile  of  an  in  ompr<  ; ill 
viscous  fluid,  as  it  flows  between  two  parallel  plates  gradually  changes, 
as  it  develops  fully  from  the  initial  profile  to  the  classical  parabolic 
profile,  seen  in  Poisouille  flow.  But,  regardless  of  the  shape  of  the 
duct’s  cross  section  no  exact  solution  is  had  for  the  profile  because  of 
the  nonlinear  advection  terms  m the  equations  of  motion  of  the  problem. 
Certain  approximate  analytic  solutions  have  been  tried  by  various  inve  t i- 
gators  using  ••ither  i series  approximation  or  a division  of  thi  ontr-mo 
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r< • |ioi.  into  two  zones.  Although  ttiere  is  more  recent  work  the:  method 

49 

b rmed  by  Schlichtmg  is  basic  to  subsequent  studies  in  that  at  the 
zone  near  the  channel  entrance,  the  boundary  layer  equations  are  applied 
and  an  approximate  solution  is  obtained,  in  terms  of  a perturbation 
serie;  . In  the  zone  where  the  profile  is  almost  fully  developed,  the 
solution  is  obtained  as  a perturbation  of  the  fully  developed  velocity 
profile.  These  two  solutions  are  matched  at  a suitable  point  to  form 
a complete  solution.  Since  this  problem  of  evaluating  the  true  velocity 
profile  at  the  duct  entrance  region  is  a difficult  one  which  requires 
more  study,  at  this  point  of  the  continuing  investigation,  fully  developed 
flow  in  the  form  of  a parabolic  velocity  profile  at  the  entrance  region, 
was  assumed. 

Theoretically,  for  a conduit  with  an  abrupt  expansion,  uniform  flow 
will  only  be  obtained  at  infinity  in  either  direction  from  the  expansion. 
Our  situation  is  assumed  to  be  a flow  into  an  expansion  region,  especially 
when  viewing  the  blood  flow  as  coming  through  the  heart  valve  seat  into 
the  sinus  of  Valsalvae  region,  upstream  of  the  obstacle.  In  turn,  the 
velocity  profile  of  uniform  laminar  flow  is  known  to  be  parabolic.  Ob- 
viously, for  a computational  solution,  the  length  of  the  conduit  must 
be  as  short  as  is  compatible  with  desired  accuracy  and  limitation  of 
computer  storage  and  time.  For  a low  Reynolds  number  the  velocity  pro- 
file approaches  a parabolic  form  after  a relatively  short  length  from 

the  entrance  in  a uniform  conduit,  according  to  the  classical  approaches 

50  51 

of  Rrandtl  and  Tietjens  and  Schlichting  . For  our  range  of  Rey- 
nolds numbers  it  was  therefore  assumed  that  a parabolic  velocity  dis- 
tribution at  a proper  finite  distance  upstream,  as  well  as  downstream 
from  the  entrance  region,  was  reasonable  until  more  detailed  invest i- 

gat  ions  could  l>  formed.  liven  for  a distance  of  3/4  b (D  being  the 

o o 
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diameter  of  the  inlet)  upstream,  Hung  has  confirmed  that  uniform  flow 
at  the  entry  region  position  would  not  be  disturbed  noticeably  by  any 
nonuniformity  on  the  downstream  side.  (Note  that  this  assumption  is 
known  to  be  untrue  for  Re  < 50  and  questionable  between  50  •'Re  -^250.) 

The  velocity  and  stream  function  profiles  of  the  fluid  leaving  the 
channel  is,  in  general,  unknown,  and  moreover,  it  is  desirable  not  to 
detail  the  profiles  explicitly.  Small  changes  in  the  exit  velocity  [Pro- 
file would  produce  large  effects  on  the  interior  stream  function  and 
vorticities.  In  order  to  allow  this  profile  to  develop  in  a natural 
fashion  rather  than  to  force  it  to  predetermined  conditions,  the  vorticity 
and  stream  function  at  the  exit  were  made  equal  to  their  values  on  the 
adjoining  interior  column  of  grid  points.  In  a rapidly  changing  velocity 
field  a time  extrapolation  was  also  included,  based  upon  changes  in  the 
vorticity  and  stream  functions  during  the  preceding  two  time  steps. 

In  some  calculations,  some  perturbation  effects  were  placed  upstream 
of  the  obstacle  in  both  a symmetric  and  asymmetric  fashion,  however,  such 
attempts  were  exploratory  ones. 

C.  Display  of  Computer  Graphics  Results  - Obstacle  and  Mesh  Lines 
Coincidental . 

Kvery  parameter  in  the  computer  program  was  kept  constant  except 
the  obstacle  shape  in  a particular  experiment,  in  order  to  determine  the 
effect  of  the  obstacle  shape  on  the  development  of  the  flow  behind  the 
obstacle . 

number  of  mesh  rows:  26 

number  of  mesh  columns:  61 

Reynolds  number:  1000  (unless  otherwise  stated) 
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error  tolerance: 

vorticity  function:  0.001  R 

e 

stream  function:  0.00003  R 

e 

or 

0.0005  R 
e 

The  series  of  computer  program  subroutines  formed  at  this  junction 
of  the  investigation  allowed  two  interesting  evaluations.  For  each 
obstacle  of  interest,  computer  graphics  results  were  developed  for  the 
stream  function,  vorti  lty  function,  horizontal  and  vertical  velocity 
components,  total  velocity  and  velocity  vector.  In  turn,  all  such  plots, 
except  the  velocity  vector,  consisted  of  two  parts;  a solutior  space 
surface  plot  and  its  isometric  projection.  An  additional  program  sub- 
routine allowed  the  presentation,  by  choice,  of  a particular  time  frame 
of  any  of  the  above  desired  functional  plots,  as  formed  with  the  aid 
of  a Gerber  plotter.  Figure  31  shows  an  example  of  a console  display 
result;  in  this  case  a stream  function. 

For  the  situation  of  the  stream  function,  such  a plot  can  be  veri- 
fied by  a laboratory  flow  visuali zat ion  technique  such  as  the  injection 
of  a dye  or  radio-active  particles  into  the  fluid.  The  stream  function 
isometric  plot  shows  the  relative  magnitude  of  the  stream  lines,  thus 
forming  a semi -quantitati ve  evaluation.  The  vorticity  function  plot 
indicates  the  magnitude  of  the  angular  momentum  of  a fluid  particle  in 
the  solution  space.  This  plot  cannot  be  obtained  by  laboratory  experi- 
ment. In  the  case  of  the  velocity  plots  the  solution  space  surface  plots 
are  equivalent  to  those  obtained  experimentally  in  the  laboratory. 

These  are  gained  by  traversing  in  a column  or  row  direction  the  space 
with  a hot  wire  anemometer  probe  so  as  to  obtain  the  magnitude  of  the 


velocity  or  desired  velocity  component.  The  contour  plot  forms  the 
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trace  of  the  equivalue  line  or  isovelocity  line  of  flow.  Such  a plot 
would  be  required  to  move  through  the  solution  space  in  such  a manner 
i!  ' ; ■ reading  held  it  a constant  value.  rhe  v<  ■ Lty  v<  • r 

plot  indicates  both  the  direction  and  relative  magnitude  of  the  flow. 
This  is  analogous  to  the  rather  crude  laboratory  method  of  visualizing 
■ II  ■ liri  :t  ion  by  inserting  "streamers"  into  the  flow  stream. 

Figure  32  is  a comparison  of  the  u-velocity  profile  for  a two 
dimensional  channel  flow  situation  as  formed  analytically  in  reference 
3.  The  significance  is  that  only  the  profile  of  the  flow  entering  the 
hannel  was  specified.  The  computer  program  calculated  its  own  • xit 
pi  file.  This  is  important  since,  in  general,  an  investigator  knows 
only  the  profile  of  the  fluid  that  enters  the  field  of  interest. 


Fi  i.  31  : A Typical  Result  For  Flow  About  An  Obstacle  As  Seen  on  the 

tDI  Scope. 
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Figure  33  shows  the  stability  values  of  the  numerical  routine  Po 
certain  stream  and  vorticity  functions  and  largest  error  positions,  as 
originally  displayed  on  the  IDI  display  scope  for  changing  iteration 
values.  Figure  33,  however,  is  a hard  copy  of  one  time  frame  within  a 
changing  time  series  seen  on  the  display  screen.  The  top  left  plot  in 
this  particular  diagrn  shows  the  anount  of  error  versus  amount  of  re- 
quired iterations  performed  by  the  computer  for  forming  a particul  ar  stra;' 
function  exhibit.  During  each  iteration  the  computer  monitored  the  al- 
lowable amount  of  error  during  the  computational  process.  When  the 
analysis  searcher  (computer)  arrived  at  the  required  error  tolerance, 
it  initiated  evaluation  of  the  stream  function  for  the  next  time  increment. 
The  middle  top  plot  shows  the  convergence  of  the  numerical  routine  to  a 
desired  value,  accomplished  by  monitoring  the  stream  function  error 
term.  A like  error  tolerance  evaluation  is  shown  for  the  vorticity 
function.  The  bottom  view  shows  positions  of  the  largest  error  terms 
for  the  vorticity  function  in  an  idealized  orthogonal  renal  section  at 
a particular  instant  of  analysis.  These  positions  are  where  the  vorticity 
function  changed  most  rapidly  during  each  iteration,  dependent  upon  para- 
metric value  listed  below  the  plot. 

The  rather  large  series  of  computational  results  on  tin  succeeding 
pages  depict  various  end  products  as  based  upon  the  previous  ly  noted  num- 
erical analysis  and  resulting  algorithms.  It  is  to  be  noted  that  all 
results  at  this  point  are  based  upon  the  requirement  that  the  boundaries 
of  the  obstacles  and  the  boundaries  of  the  computing  mesh  cello  are  coin- 
cidental. Close  observation particularly  for  the  oval  and  ciieular 
; ■ icles,  thi  >W! i that  th<  ■■  ol  je  t . > 1 boundarii  e >r<  formed  by 

lines;  at  the  most,  any  indivi  iuul  line  crossing  a mesh  cell  was  t li agon,, 
>n<  , allowed  by  th<  algorithm.  Curved  boundaries  are  disc  us  ••  1 in  tic 
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Flcj.  31:  . Lability  Valuta  of  the  Computers  Numerical  Routine  for 
Certain  Stream  and  Vorticity  Functions.  Positions  of 
Largest  error  Terms  for  a Vorticity  Function  with  a Par- 
ticular Reynolds  Number  Value  is  also  Shown. 
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Each  set  of  figures  for  the  rectangular  (figures  34-48),  - 

hemisphere  (figures  49-62),  circular  (figures  63-78)  and  oval  shaped 
(figures  79-93)  obstacles  are  required  for  the  prosthetic  heart  valve 
section  of  the  investigation.  Particular  stream  function  franv  s ar 
first  shown  to  point  out  the  capability  of  the  graphics  display  for 
allowing  study  of  a continuously  changing  flow  field  (remembering  that 
the  present  requirement  is  blood  movement  downstream  of  a heart  valv 
occluder  or  about  an  atherosclerotic  lesion).  Results  shown  are  hard 
copy  equivalences  to  the  console  display.  Subsequently,  examples  of  the 
u and  v and  total  velocity  plots,  as  well  as  the  velocity  vector  plot  are 
shown  for  such  obstacle  shapes.  The  velocity  vector  plot  allows  the 
added  advantage  that  if  the  velocity  vector  arrows  are  not  consistent 
with  their  neighbors’  paths  (including  circular  motion  for  the  vorti  ity 
positions)  it  is  a warning  that  there  may  be  a computational  error  at 
that  position.  Arbitrary  Reynolds  number  values  are  shown  in  each  figure 

Another  set  of  figures  involve  a hexagonal  obstacle  (figures  34-104) 
and  its  functional  plots  for  two  different  time  frames. 

As  an  introduction  to  the  possibility  of  employment  of  computer 
graphics  for  studying  valve  leaflet  motion,  one  simplified  leaflet  form 
in  a steady  flow  field  is  shown  in  figures  105-112  for  two  different 
Reynolds  numbers  and  time  values.  Figures  113-118  show  evaluation  •.  <>s- 
sibilities  for  flow  from  the  atrium  to  the  ventricle  section  of  the  heart 
in  a very  rudimentary  shape  arid  fashion,  for  the  present.  In  all  oases 
oru  views  rigid  wall  and  steady  flow  conditions.  This  point  will  be 
|j  > ussed  in  a Latei  section.  In  turn , subroutines  t < > t >rm  figuri  11  I 
ip  introduction:  to  a two  leaflet  natural  valve  study  or  ot.het  t tot  bo 
that  require  similar  shaped  sections. 
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Figures  121-1 31  again  show  other  usages  for  computer  graphics.  In 

tins  < xample,  one  time  frame  plots  are  shown  for  two  obstacles  in  a flow 

* ★ 

field.  One  Last  example  in  this  section  is  the  stream  function  pi  >t  of 
figure  132  which  stresses  the  capability  of  the  procedure  as  it  could 
pertain  to  a computer  graphics  analysis  of  an  arbitrary  aircraft  wing 
shape  placed  at  any  flow  angle. 


**  This  study  was  evolved  by  Mr.  Charles  Brauer  for  a particular  fluidic 
; • tudy  . required  by  thi  Fluidonics  Corp.  il  ill  Lak<  :ity 
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Fig.  34:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a Rectangular 

Obstacle 
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Figure  35:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a Rectangular 

Obstacle 
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Fig.  36:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  37:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  38:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  39:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  40:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  41:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectnagular  Obstacle. 
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Fig  42:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Rectangular  Obstacle. 
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Fig.  43:  Stream  Function  Plot  at  a Particular  Time  Trame  for  a 

Rectfna^ular  Obstacle. 
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'eciit.  Function  Plot  at  a Particular  Time  Frame  for  a Rectangular 
• ey  no 1 d . . . . . 
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Fiq.  45:  The  Horizontal  Velocity  Function  Plot  for  a Rectanqul... 

Obstacle  (Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  49:  Stream  Function  Plot  at  a Particular  Time  Frame  for 

Semi-Hemisphere  Obstacle. 
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Fig.  50  : Stream  Function  Plot  at  a Patticular  Time  Frame  for 

Semi-Hemisphere  Obstacle. 
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Fig.  51:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi -Hemisphere  Obstacle. 
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Fig.  52:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  53:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  54:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  55:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  56:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  57:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  58:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Semi-Hemisphere  Obstacle. 
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Fig.  60:  The  Vertical  Velocity  Function  Plot  for  a Semi -Hemisphere 
Obstacle  (Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  61:  The  Total  Velocity  Function  Plot  for  a Semi-Hemisphere  obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  62:  The  Velocity  Vector  Plot  for  a Semi-Hemisphere  Obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  63:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 

I 


I* 


H K:00. 

j l-j  'i  ’ 


Fig.  64:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  65:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  66:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  67:  Stream  Function  Plot  at  a Particular  Time  Frame  for 

Spherical  Obstacle. 


I J I fU  I 


Fig.  68:  Stream  Function  Plot  at  a Particular  Time  Frame  for 

Spherical  Obstacle. 
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Fig.  69:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  70:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  71:  Stream  Function  Plot  at  a Particular  Frame  for  a 

Spherical  Obstacle. 
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Fig.  72:  Stream  Function  Plot  at  a Particular  Time  Frame  for  a 

Spherical  Obstacle. 
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Fig.  74:  The  Vertical  Velocity  Function  Plot  for  a Spherical  Obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  75:  The  Total  Velocity  Function  Plot  for  a Spherical 

Obstacle  (Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  76:  The  Velocity 
(Surface  and 


Vector  Plot  for  a Spherical  Obstacle 
Isometric  Views  at  One  Time  Frame) 


Fia.  77:  An  Experimental  View  of  Flow  Streamlines 

About  a Spherical  Obstacle  (From  Fromm40 
Flow  is  From  Right  to  Left. 
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Fig.  79:  Stream  Function  Plot  at  a Particular  Time  Frame  for  An 

Oval  Obstacle. 
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Fiu.  HO;  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

oval  Obstacle. 
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Fig.  81:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  82:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  83:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  84:  Stream  Firction  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  85:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  86:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  87:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 


Fig.  88:  Stream  Function  Plot  at  a Particular  Time  Frame  for  an 

Oval  Obstacle. 
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Fig.  89:  The  Horizontal  Velocity  Function  Plot  for  an  Oval  Obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  90:  The  Vertical  Velocity  Function  Plot  for  an  oval  Obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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Fig.  91:  The  Total  Velocity  Function  Plot  for  an  Oval  Obstacle 
(Surface  and  Isometric  Views  at  One  Time  Frame) 
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A Stream  Function  Plot  at  a Particular  Time  Frame  and 
High  Reynolds  Number  for  an  Oval  Obstacle  (Frames  41  and  44) 
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Fig.  94:  The  Stream  Function  Plot  for  an  Octagonal  Obstacle 
(Surface  and  Isometric  Views) 
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Fig.  95:  The  Stream  Function  Plot 
for  an  Octagonal  Obstacle 
(Surface  and  Isometric 
Views  at  a Later  Time. 
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Fig.  97:  The  Horizontal  Velocity  Function  Plot  for  an  Octagonal 
Obstacle  (Surface  and  Isometric  Views  at  a Later  Time) 
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Fig.  98:  The  Vertical  Velocity  Function  Plot  for  an  Octagonal  Obstacle 
(Surface  and  Isometric  Views) 
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Fig.  99:  The  Vertical  Velocity  Function  Plot  for  an  Octagonal  Obstacle 
(Surface  and  Isometric  Views  at  a Later  Time) 


Ill 


t o r.uL  vnincirr 
nt  man. 
muni-;  m.  m. 

Fig.  100s  The  Total  Velocity  Function  Plot  for  an  Octagonal 
(Surface  and  Isometric  Views) 
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Fig.  101:  The  Total  Velocity  Function  Plot  for  an  Octagonal 
(Surface  and  Isometric  Views  at  a Later  Time) 
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Fig.  102:  The  Vorticity  Function  Plot  for  an  Octagonal  Obstacle 
(Surface  and  Isometric  Views) 
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Fig.  103:  The  Vorticity  Function  Plot  for  an  Octagonal  Obstacle 
(Surface  and  Isometric  Views  at  a Later  Time) 
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104:  The  Velocity  Vector  Plot  for  an  Octagonal  obstacle 

(Surface  and  Isometric  Views  at  a Particular  Time  Frame) 
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Fig.  105:  Stream  Function  Plot  for  Angled  One-Cusp  Rigid  Heart 
Valve  Leaflet  (One  Time  Frame) 


Fig.  106:  Stream  Function  Plot  for  Angled  One-Cusp  Rigid  Heart  Valve 
Leaflet  (One  Time  Frame  and  at  a Higher  Reynolds  Number) 
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Fig.  107:  Vorticity  Function  Plot  for  Angled  One-Cusp  Rigid  Heart 
Valve  Leaflet  (One  Time  Frame) 
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Fig.  108:  Horizontal  Velocity  Function  Plot  for  Angled  One-Cusp  Rigid 
Heart  Valve  Leaflet  (One  Time  Frame  and  at  a Higher  Reynolds 
Number) 
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Fig.  109:  Vertical  Velocity  Function  Plot  for  Angled  One-Cusp  Rigid 
Heart  Valve  Leaflet  (One  Time  Frame) 
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Fig.  110:  Total  Velocity  Function  Plot  for  Angled  One-Cusp  Rigid 

Heart  Valve  Leaflet  (One  Time  Frame  and  at  a Higher  Reynolds 
Number) 
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Fig.  Ill:  Velocity  Vector  Function  Plot  for  Angled  One-Cusp  Rigid  Heart 
Valve  Leaflet  (One  Time  Frame) 
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Fig.  112:  Velocity  Vector  Plot  for  Angled  One-Cusp  Rigid  Heart 

Valve  Leaflet  (One  Time  Frame  and  Higher  Reynolds  Number) 


Fig. 


113:  Stream  Function  Plot  for  Idealized  Flow  into  Simplified 
Atrium  Section  of  the  Heart  (One  Time  Frame) 
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Fig.  114:  Horizontal  Velocity  Function  Plot  for  Idealized  Flow  into 
Simplified  Atrium  Section  of  the  Heart  (One  Time  Frame) 
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Fig.  116:  Total  Velocity  Function  Plot  for  Idealized  Flow  into  Simplified 
Atrium  Section  of  the  Heart  (One  Time  Frame) 
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Fig.  117:  Vorticity  Function  Plot  for  Idealized  Flow  into  Simplified 
Atrium  Section  of  the  Heart  (One  Time  Frame) 
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Fig.  118:  Velocity  Vector  Plot  for  Idealized  Flow  into  Simplified 
Atrium  Section  of  the  Heart  (One  Time  Frame) 
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Fig.  119:  Stream  Function  Plot  for  Two  Arbitrary  Shaped  Flow  Areas 
(Amenable  to  Particular  Applications  such  as  Two  Cusped 
Leaflet) 
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Fig.  120:  Horizontal  Velocity  Function  Plot  for  Two  Arbitrary  Shaped 
Flow  Areas  (One  Time  Frame) 


V VELOCITY 
Rr  U(lO. 

FRMMC  NO.  10. 


Fig.  121:  Vertical  Velocity  Function  Plot  for  Two  Arbitrary  Shaped  Flow 
Areas  (One  Time  Frame) 
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Fig.  122:  Total  Velocity  Function  Plot  for  Two  Arbitrary  Shaped  Flow 
Areas  (One  Time  Frame) 
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Fig.  123:  Vorticity  Function  Plot  for  Two  Arbitrary  Shaped  Flow  Areas 
(One  Time  Frame) 
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Fig.  124:  Velocity  Vector  Plot  for  Two  Arbitrary  Shaped  Flow  Areas 
(One  Time  Frame) 
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Fig.  125:  Stream  Function  Plot  for  Two  Obstacles  in  a Flow  Field 
(One  Time  Frame) 
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Fig.  126:  Photograph  fro.:.  IDI  Console  Screen  Showing  Stream  Function 
Plotted  for  Detailed  Studv  of  Internal  Core  of  Flow. 
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Fig.  127:  Horizontal  Velocity  Function  Plot  for  Two  Obstacles  in 
Field  (One  Tine  Frano) 
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Pig . 128:  Vertical  Velocity  Function  Plot  for  Two  Obstacles  in  a Flow 
Field  (One  Time  Frame) 
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Fig.  129:  Total  Velcoity  Function  Plot  for  Two  Obstacles  in  a Flow  Field 
(One  Time  Frame) 
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Fig . 130:  Vorticity  Function  Plot  for  Two  Obstacles  in  a Flow 
Field  (One  Time  Frame) 
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131:  Velocity  Vector  Plot  for  Two  Obstacles  in  a Flow  Field 
(One  Time  Frame) 
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Fig.  132:  Stream  Function  Plot  for  An  Arbitrarv  Aircraft  Wing  Shape 
(One  Time  Frame) 
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CHAPTER  VI 

PRESENTATION  OF  A CURVED  BOUNDARY  ALGORITHM 

The  main  objection  to  the  previous  attempts  at  a numerical  solution 
to  the  general  set  of  the  Navier-Stokes  equations  is  that  they  are  not 
general.  Previous  investigators,  as  well  as  the  present  ones  up  to  this 
point,  assumed  that  the  boundary  of  the  finite  difference  mesh  constituted 
either  a fluid  boundary  or  a solid  boundary  that  coincided  with  the  com- 
putational mesh  lines.  Many  problems  are  ignored  and  the  computer  pro- 
gram is  greatly  simplified  if  one  tries  to  fit  all  physical  boundary  con- 
figurations to  the  lines  of  a square  finite  difference  mesh,  however, 
the  usefulness  of  such  a computer  program  is  doubtful.  In  general, 
most  flow,  design,  and/or  obstacle  problems  require  an  irregular  boundary 
configuration.  The  algorithm  presented  in  this  section  is  an  attempt 
to  form  the  required  numerical  routine  for  allowing  curved  boundaries 
into  the  desired  two  dimensional  flow  problem.  The  finite  difference 

equations  for  forming  the  curved  boundary  problem  were  thus  formed. 

54 

Specialized  finite  difference  concepts  are  discussed  by  Forsythe  , 

• i j .55  56 

salvadon  , and  Todd 

The  computational  mesh  with  stencils  of  unequal  arms,  as  discussed 
in  section  B2,  was  employed,  as  well  as  the  finite  difference  forms  of 
section  Bi.  As  before,  the  following  definitions  were  set; 

a.  Let  the  superscript  k denote  the  value  of  the  variables 
at  the  kth  iteration. 

b.  Lot  the  subscripts  A,B,C  denote  the  boundary  points, 
adjacent  interior  points,  and  remaining  interior  points, 
respectively,  of  the  node  points  on  the  finite  difference 


IUI-.'D 


k 1 1 

mesh  grill,  as  seen  in  figure  13  3.  For  « xample  t| ^ , 


H Adjacent  inte  rior  point  ( t>  pe  0) 
A Interior  point  (type  0) 

© Node  boundary  point  (type  1) 
Fxterior  node  point  (type  2) 


Computational  Mesh  with  Curved  Boundary  Involved. 
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denotes  the  value  of  the  stream  function  evaluated  at 
the  k + 1th  iteration  on  the  point  (i,j)  where  (i,j) 
corresponds  to  a boundary  point . 

The  algorithm  follows  the  previous  one  up  to  a certain  point.  For 
the  sake  of  continuity,  the  totality  of  steps  are  shown. 

1.  The  values  il  , w were  estimated  by  extrapolation  from 

13  j C B | C 

past  values. 

2.  The  improved  values  of  were  calculated  from 


m = v u)  - b uj  + ib  uj 

t V x x y 


(in  finite  difference  form,  with  the  right  hand  side  averaged 
between  level  k and  k + 1) . 

Again,  the  method  of  calculation  was  the  Peaceman-Rachford 
alternating  direction  method.  Smoothing  was  neeess'ary  for 
higher  Reynolds  numbers. 

3.  The  stream  function  values  on  the  adjacent  interior  points 
were  improved  by  the  following  formulae; 

*B2  * *<W  W * hV VVW 

* -yy  <vvw  + 0lh3>  1591 

Since  = -2ui  - ip  then 

yy  xx 


b = — 2oj  — — ( ib  —2ib  + ib  ) 

yy  „ b ,2  o V 

11  B h 


But  is  a boundary  point  and  was  not  known,  therefore 
extrapolation  from  known  interior  points  was  required, 


tij  — <*) 

H B2  B * 
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and  the  following  improvement  formula  was  required; 


k+1  k u k h 

- il>  + h ti;  + — 

^B2  vyB  2 


„ k „ k 1 
- 4uj  + 2uj  - — 
B2  B3  u2 
h 


_ k k ) 

2 t/i  + ili 

B C 


(61) 


Since  the  boundary  conditions  \p*  = ip^  = 0 were  necessarily 
satisfied  everywhere  on  a fixed  solid  boundary,  then  if  the 
boundary  is  parallel  to  the  mesh, 


k 

if;  _ = 0 

yB 

and  equation  (61)  reduced  to 
, 2 


k+1  k 

\1  „ = ili  + — 

u32  VB  2 


- 4d>. 


B2 


2uj. 


B 3 


k „ k k 

ili  - 2 ip  + ip 

A B 


(62) 


If  the  boundary  is  curved,  however,  as  depicted  in  figure  133  (line  C ) 

then  one  must  make  some  additional  calculations.  Consider  the  mesh  node 

point  (i,j)  in  figure  133.  The  net  boundary  is  shown  by  the  dashed  line 

. At  each  node  of  C^,  a perpendicular  was  drawn  through  and  extended 

until  it  crossed  a link  or  node  of  the  net  in  the  interior.  For  example, 

the  perpendicular  from  B meets  C at  F,  and  meets  the  link  from  B to  C_, 

B 2 2 

involving  the  angle  u. 

A value  for  ip  was  determined  by  linear  interpolation 
between  ip^^  and  A "normal"  derivative  for  the 

linear  interpolation  was  then  obtained  for  the  net 
function  ip  , by 


h/cos  a 


(6  3) 


14‘3 


Equating  (63)  to  the  normal  derivative  at  F,  one  obtained 

the  desired  equation  for  . For  example,  expressing  j/ 

in  terms  of  6 and  <p  by 
B z Lz 


B2»E  = h tan  a,  E«  C2  = h(  1-tan  a) 


(64  a,b) 


*B2(cos  a “ sin  + '^c2  sin  a ~ J'B  cos  a = 


hi(-  =0 

" F 


hence , 


ip  = — ip  _ (cos  a -sin  a)  + if  sin  a)  (66) 
B cos  a B 2 C2 


The  equation  for  the  adjacent  interior  node  B2  that  was 
analogous  to  equation  (61)  was  then 


k+1  k h 

lb  = lb  + — — 

tB2  rB  2 


, k , k l 

%2  - 2U)B3  " ~2 
h 


, k _ k kv 

( ill  - 2ip  + <p  ) 

1 B 


k k k 

where  , and  are  pseudo-boundary  points  that 

were  obtained  by  extrapolation  from  the  known  values,  and 

m- 

k k k 

where  if  , ip  , and  ip  all  satisfied  equation  (66)  . 

A B C 

4.  The  improved  values  of  ui  were  calculated  by 

B 

NEW  (i)R  = SF  » (OLD  t»>B ) + (1-SF)  T 
where  T = -l/2\’’ip\  and  where  SF  is  a smoothing  factor. 


5.  Test  for  convergence;  if  not  convergent  one  returned  to  step  1. 


\ 
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A.  Discussion  for  Atherosclerosis  Study 

Examples  of  computer  graphic  flow  plots  follow.  All  are  equivalent 
to  the  original  "experiments"  as  performed  on  the  display  scope,  however, 
the'  alternate  program  sub-routine  for  employing  the  Gerber  plotter  was 
used  to  allow  the  particular  one  time-frame  results  shown.  Parametric 
flow  solutions  are  first  shown  for  particular  vascular  situations  in  the 
renal  section  and  in  the  abdominal  aortic  bifurcation  section.  This 
group  is  then  followed  by  graphic  results  formed  for  certain  artificial 
heart  valve  situations. 

In  most  evaluations  the  plots  for  the  stream  function,  horizontal 
velocity  component,  vertical  velocity  component,  total  velocity  function, 
vorticity  function  and  velocity  vector  are  shown.  All  have  their  surface 
and  isometric  values  formed.  Although  details  of  concern  to  each  plot 
have  not  been  pursued,  they  are  available  for  future  study.  One  can, 
however,  by  qualitative  viewing,  infer  some  interesting  points.  For 
example,  figures  135  and  137  show  the  typical  laminar  flow  front,  to  be 
expected  at  the  Reynolds  number  employed  in  those  results.  In  turn,  the 
isometric  plots  of  figure  135  and  136  combine  to  form  the  isometric  plot 
in  figure  137.  Flow  is  from  left  to  right.  In  some  instances  the  hori- 
zontal velocity  plot  appears  to  be  equivalent  to  the  total  velocity  plot 
in  isometric  view.  Qualitatively  this  is  equal  to  saying  that  the  vertical 
velocity  component  was  extremely  small  in  value.  Figure  138  shows  heavy 
concentrations  of  vortices  (forerunners  of  turbulence)  rapidly  appearing 
at  the  entrances  to  the  smaller  arterial  channels,  with  interacting 
vortices  after  a right  angle  turn  from  the  main  channel  flow.  The 
equivalent  isometric  picture  snows  the  sudden  rises  in  vortex  magnitudes 
it  these  clustered  positions.  The  isometric  plot  allows  quantitative 
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evaluation  since  the  vortex  peaks  arise  from  a flat  background  that 
constitutes  a baseline  value. 

Note  that  fiqure  134  through  139  did  not  require  the  curved  boundary 
algorithm.  Figure  140,  as  an  example  of  the  first  use  of  the  curved 
boundary  algorithm,  does  show  a heavy  vortex  concentration  within  the 
upper  side  artery  leading  to  one  kidney  from  the  mesentery  artery  theo- 
retical equivalent.  Depending  upon  the  particular  time  frame  employed, 
if  viewed  in  the  light  cf  a Von  Karman  vortex  street  likeness,  that 
particular  vortex  could  have  appeared  in  the  other  side  artery. 

Figures  141  through  145  depict  computer  graphic  results  as  applied 
to  the  more  natural  situation  of  non-orthogonal  side  channels  emanating 
from  the  main  artery  in  the  renal  section.  The  important  point  is  that 
the  side  arteries  are  programmed  to  angle  at  any  desired  inclination, 
just  as  clinical  results  show  for  various  individuals.  The  upper  side 
artery  was  arbitrarily  placed  facing  upstream  and  the  lower  one  was 
faced  downstream  to  the  general  flow. 

Figure  141  bears  out  physical  concepts  in  that  the  streamlines  are 
shown  closer  in  proximity  upstream  of  the  side  arteries,  therefore, 
showing  lessening  of  flow  energy  downstream  as  flow  is  drained  off  in 
the  side  channels.  Figure  142  shows  the  typical  parabolic  type  front 
for  laminar  flow,  changing  to  a typical  turbulent  front  in  the  main 
flow  region  between  the  side  arteries.  Figure  144  shows  a proper 
physical  amalgamation  of  figures  142  and  143,  the  usual  horizontal 
and  vertical  velocity  components  together  forming  a total  velocity 
value  in  the  isometric  view.  Figure  145  shows  the  vorticity  function 
which  then  became  tin;  subject  of  d<  tailed  sequential  study  within  the 
particular  portion  of  the  backward  facing  side  artery  region,  resulting 


in  figure  140. 


Figure  145's  isometric  plot  particularly  shows  the  extremely 
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hiqh  vortex  peaks  in  the  upper  channel,  which  was  to  be  expected  as  the 
blood  flow  radically  changed  direction.  Realizing  this,  the  Reynolds 
number  was  continually  changed  in  the  montage  components  of  figure  146 
to  allow  a transition  from  a slow  initial  rate  of  flow  to  slightly  above 
an  average  blood  flow  velocity  value  in  that  particular  area.  One  then 
notes  the  formation  of  a vortex  in  such  a velocity  regime  and  if  the  ob- 
server views  the  small  triangular  shaped  vortex  on  the  anterior  wall, 
the  breakup  and  reformation  of  stream  function  lines  in  amplifying  this 
vortex  are  clearly  noticed. 

Figure  117  depicts  the  various  plots  as  formed  for  an  idealized 
atherosclerotic  lesion  formed  on  a vascular  wall.  Present  work,  involving 
stress  analysis  at  the  vascular  wall,  but  not  reported  here,  will  certainly 
incorporate  the  program  subroutines  that  formed  this  particular  picture. 

This  study  area  is  to  be  extended  and  one  will  attempt  to  understand  the 
growth  and  interplay  of  flow  parameters  as  the  plaque  grows. 

Figures  148  through  153  show  the  stream  function,  horizontal  velocity 
component,  vertical  velocity  component,  total  velocity,  vorticity  function 
and  velocity  vector  for  a particular  Reynolds  number  value  at  the  abdominal 
aortic  bifurcation  region.  Figure  152  may  be  singled  out  as  particularly 
interesting  and  for  comparison  to  figure  154. 

Note  that  figure  154  is  a diagrammatic  representation  of  sites  of 

23 

predilection  for  the  occurrence  of  atherosclerosis.  Wesolowski  and 
other  vascular  suigeons,  using  hydraulic  models  of  various  configurations, 
and  from  observations  at  operations  and  by  arteriography,  have  noted  that 
turbulence  occurs  at  each  site  of  such  predilection.  Comparison  then,  of 
figure  152  with  figure  154  shows  that  the  theoretical  vorticity  function 
plot  (vorticity,  the  forerunner  of  turbulence),  is  equivalent  for  site 
positioning  with  observed  sites  of  turbulence  and  eventually  atherosclerosis 


site  prediction. 
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Figures  155  through  158,  and  160  and  161,  view  the  computer  graphic 

results,  at  a particular  instant  of  time  and  for  an  arbitrary  Reynolds 

number,  for  blood  flow  within  the  renal  and  abdominal  aortic  bifurcation 

sections  together.  The  physical  flow  values,  consistent  with  those  of  an 

average,  healthy,  young  male,  will  be  introduced  into  such  a computer 

model  in  the  near  future.  If  one  compares  figure  160  with  figure  152,  it 

appears  that  the  vortices  at  the  bifurcation  are  shifted  downstream  in 

figure  160  and  added  lengths  of  vortex  motion  are  noted.  Such  results 

are  hypothesized  to  be  due  to  the  role  played  by  the  renal  section  as 

5 7 58 

an  upstream  addition.  This  concept  is  discussed  by  Greenfield  ' 

Figure  160,  when  compared  to  figure  161,  shows  that  orthogonal  side 
arteries  at  the  renal  section  do  have  somewhat  different  vortex  motion 
than  askew  arteries.  In  addition,  however,  the  angled  renal  arteries 
appear  to  produce  heavier  vortex  motion  within  the  bifurcated  channels 
downstream. 
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B.  Curved  Boundary  Computer  Graphic  Displays. 


(Vascular) 


Fig.  134:  The  Stream  Function  Plot  for  Idealized  Renal  Artery  Section 
(Orthogonal  Branching  Arteries;  One  Time  Frame) 


135:  Horizontal  Velocity  Function  Plot  for  Idealized  Renal  Artery 
Section  •(Orthogonal  Branching  Arteries;  One  Time  Frame) 
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136:  Vertical  Velocity  Function  Plot  for  Idealized  Renal  Artery 
Section  (Orthogonal  Branching  Arteries;  One  Time  Frame) 
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137:  Total  Velocity  Function  Plot  for  Idealized  Renal  Artery 
Section  (Orthogonal  •Branching  Arteries;  One  Time  Frame) 
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Fiq.  139:  Velocity  Vector  Riot  for  Idealized  Renal  Artery  Section 
(Orthogonal  Branching  Arteries;  One  Time  Frame) 
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Fig.  141:  Stream  Function  Plot  for  Branching  Arteries  at  Renal  Section 
(One  Branch  Faces  Upstream  and  the  Other  Downstream)  at  the 
Renal  Section.  One  Time  Frame. 
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Fig.  142:  Horizontal  Velocity  Function  Plot  for  Branching  Arteries  at 

Renal  Section  (opposite  Facing  Branch  Arteries) . One  Time  Frame 


143:  Vertical  Velocity  Function  Plot  for  Branching  Artereis  at 
Renal  Section  (One  Branch  Artery  Faces  Downstream  and  One 
Faces  Upstream).  One  Time  Frame. 
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144:  Total  Velocity  Function  Plot  for  Branching  Arteries  at  Renal 
Section  (One  Branch  Artery  Faces  Upstream  and  One  Faces 
Dows tr earn ) . One  Time  Frame. 
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FRAME  NO.  10. 

145:  Vorticity  Function  Clot  for  Branching  Arteries  at  Renal  Section 

(One  Branch  Artery  Faces  Upstream  and  the  Other  Faces  Downstream) 
One  Time  Frame. 
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Reynolds  Number  = 120 


Fig . 


146:  Montage  of  Portion  of  Renal  Section  with  Side  Artery  Facing  Upstream. 
At  Different  Velocities,  Showing  Vortex  Development. 
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Fig.  147:  A Montage  That  Shows  Various  Computer  Devised  Flow 
Functions  For  an  Idealized  Atherosclerotic  Plaque. 


U VELOCJTT 
H-  150. 

FPiRrtE  NO.  15. 

Horizontal  Velocity  Fun,  t i 
Bifurcation  (One  Tine  r, 


AD-A038  694 


UNCLASSIFIED 


UTAH  UNIV  SALT  LAKE  CXTT  OEPT  OF  COMPUTER  SCIENCE  F/8  6/5 

^«AT?LI5AIi2r^r.C2MP«TE-  GRAPHICS J TWO  CONCURRENT  INVESTI8ATX— ETC  (U) 
NOV  71  H 8REENFIEL0*  R DEBRY  F30602-70-C-0300 

UTEC-CSC-71-115  U2r° 


3 OF  £/- 

AS 038694 

■ 

> 

> 

1 ^ 
1 ^ 

s*=^ 

> 

> 

<* 

c: 

a. 

% 

LJ  1 

% 

C: 

% 

c r 

C= 

H 

42'r: 

< ~ 

TOTAL  VELOCITY 
R=  150. 

FRAME  NO.  15. 

Fig.  151:  Total  Velocity  Function  Plot  for  Abdominal  Aortic  Bifurcation 
(One  Time  Frame) 


153:  Velocity  Vector  Plot  for  Abdominal  Aortic  Bifurcation 
(One  Time  Frame) 
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Fig.  154:  Arterial  "Tree"  Showing  Sites  of  Predilection  for  Atherosclerotic 
Plaques.  (Note  Dark  Sections  Along  Walls  at  Mesenteric  Section) 
From  Wesolowski^-*. 
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Fig.  157:  Vertical  Velocity  Function  Plot  for  Renal-Nominal  Aortic 
Bifurcation  Section  Combined;  One  Time  Frame. 

V VELOCITY 
R=  150. 

FRRMi:  NO.  15. 


Fig.  158:  Total  Velocity  Function  Plot  for  Renal-Abdominal  Aortic 
Bifurcation  Combination  Section;  One  Time  Frame. 
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Fig.  159:  Vorticity  Function  Plot  for  Renal-Abdominal  Aortic  Bifurcation 
Combination  Section;  Orthogonal  Branching  at  Renal  Position. 
(One  Time  Frame) . 


Fig.  160:  Vorticity  Function  Plot  for  Renal-Abdominal  Aortic  Bifurcation 
Combination  Section;  Non-Orthogonal  Branching  at  Renal  Section 
(One  Time  Frame) . 
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Fig.  161:  Velocity  Vector  Plot  for  Renal-Abdominal  Aortic  Bifurcation 
Combination  Section.  (One  Time  Frame) 
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C.  Discussion  of  Concern  to  Prosthetic  Heart  Valve  Study. 

Figures  162  through  167  present  one  time  frame  examples  of  the  com- 
puter graphic  results  as  applied  to  an  ideal  ball  type  central  occlusion 
artificial  heart  valve.  In  turn,  figures  168  through  173  display  equivalent 
plots  for  a disc  type  artificial  heart  valve. 

Remembering  that  the  computer  is  an  inert  brute  with  non-diagnostic 
capabilities  in  the  clinical  sense,  perhaps  the  advocates  of  the  disc 
valve  will  be  pleased  to  notice  what  qualitatively  appears  to  be  larger 
vortex  formations  distal  to  the  ball  in  figure  166,  than  for  the  disc 
valve,  when  viewing  the  isometric  result.  In  comparison,  ball  valve 
proponents,  upon  viewing  figure  172,  could  perhaps  smugly  point  to  heavy 
vortex  concentration  at  the  edges  of  the  discoid,  possibly  leading  to 
quicker  edge  wear  for  this  valve  type.  With  such  a thought  in  mind,  the 
ball  type  also  appears  to  show  more  evenly  distributed  vorticity  motion 
about  the  occluder. 

Figures  175  and  176  show  the  ball  and  disc  valves  in  less  than  full 
open  position  as  the  computer  studies  them  when  the  occluder  is  closer 
to  the  valve's  sewing  position. 

Again,  qualitatively,  the  isometric  views  show  higher  peaks  for  vor- 
ticity which  is  to  be  expected  when  the  shearing  effects  are  greater 
between  the  artery  wall  and  the  occluder  as  this  space  is  lessened.  As 
one  should  also  assume  intuitively,  the  "computer  experiment"  shows  heavier 
vortex  concentration  about  both  types  of  occluders  as  they  move  towards 
closure  position.  One  can  note  the  extent  and  activity  of  the  vortex 
concentration , when  the  valve  occluder  for  either  valve  type  is  not  fully 
opened,  when  the  observer  realizes  that  figures  175  and  176  were  formed 
at  the  second  time  frame  after  initiation  of  steady  flow.  In  turn,  figures 
166  and  172  were  formed  at  the  fifteenth  time  frame.  This  section  of  the 
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investigation  is  continuing,  as  time  permits,  since  closeness  of  the 

artery  wall  requires  a finer  computational  mesh  for  the  numerical  analysis. 

Most  important  is  the  comparison  between  figures  172  and  174.  Figure 

59 

174  is  a flow  pattern  as  formed  by  Wieting  while  photographing  the 
motion  of  particles  about  a Kay-Shiley  valve  in  a pulse  duplication 
system.  Wieting's  experiments  were  performed  independently  at  Baylor 
University  and  for  preliminary  comparison  of  vortex  positioning  with  the 
present  computer  result,  the  product  is  gratifying. 
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Curved  Boundary  Computer  Graphic  Displays  (Heart  Valve) . 


STRERM  FUNCTION 
ft  150. 

FRflME  NO.  15. 

Fig.  162:  Stream  Function  Plot  for  Starr-Edwards  Type  Artificial 
Heart  Valve  at  Full  Open  Position  (One  Time  Frame) 
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FRAME  NO.  15. 

Fig.  164:  Vertical  Velocity  Function  Plot  for  Starr-Edwards  Type  Artificial 
Heart  Valve  at  Full  Open  Position  (One  Time  Frame) 
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Fig.  167:  Velocity  Vector  Plot  for  Starr-Edwards  Heart  Valve  at  Full 
Open  Position  at  Full  Open  Position  (One  Time  Frame) . 
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172:  Vorticity  Function  Plot  for  a Disc-Type  Artificial  Heart  Valve 
at  Full  Open  Position  (One  Time  Frame) 
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Fig.  173:  Velocity  Vector  Plot  for  a Disc-Type  Artificial  Heart.  Valve 
at  Full  Open  position  (One  Time  Frame) 


175:  Vorticity  Function  Plot  for  a Starr-Edwards  Type  Artificial 
Heart  Valve  at  Partially  Open  Position  (One  Time  Frame) 
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176:  Vorticity  Function  Plot  for  a Disc-Type  Artificial  Heart 
Valve  at  Partially  Open  Position  (One  Time  Frame) 
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CHAPTER  VII 


CONCLUSIONS . 

It  is  obvious  that  fluid  flow  experimentation  by  computer  graphics 
is  an  extremely  versatile  tool.  Less  cost  and  time  are  required  to  study 
design  concepts  by  computer  experiment  than  by  physical  experiments.  Also, 
in  wind  tunnel  experiments  or  in  hydraulic  modeling,  modifications  of  data 
interpretation  are  usually  required  to  compensate  for  flow  velocity  uncer- 
tainties, viscosity  variation  and  fluid  compressibility.  Such  variations 
are  easily  written  into  or  ignored  in  computer  programs.  Detailed  com- 
puter results  form  working  approximations  to  intractable  mathematical 
concepts . 

Although  the  primary  function  of  this  presentation  has  been  accom- 
plished, namely,  the  viewing  of  a powerful  new  approach  with  subsequent 
important  applications,  it  is  realized  that  a great  deal  of  effort  still 
remains  for  true  duplication  of  natural  blood  flow  in  vascular  channels. 
This  point  has  been  discussed  by  (Jreenfield  ' (>“ . The  fluidal 

system  involves  not  a steady  motion  but  a pulsatile  motion.  The  channel 
walls  of  the  arterial  tree  are  not  fixed,  solid  walls,  but  rather,  they 
display  elastomeric  or  viscoelastic  properties.  Such  properties  will 
be  seen  in  future  models  formed  for  graphic  display.  Other  effects,  such 
as  tapering  of  vascular  channels  and  the  presence  of  reflective  waves 
emanating  from  the  pulsatile  motion,  will  serve  as  future  inputs  to  the 
overall  program.  All  these  effects  will  serve  as  the  bases  for  continuing 
reports  in  this  series. 

The  computer  experiment  approach  could  replace  present  flow  pattern 
studies  formed  by  dye  observations  since  observations  cease  when  dye  and 
other  flow  components  become  mixed,  obscuring  the  flow  details.  At  the 
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Least,  t!it  computer  experiment  would  complement  the  present  tool  of 
streaming  birefringence.  Numerically,  in  the  computer  experiment,  addi- 
tional perturbation  effects  would  possibly  be  injected  into  the  flow 
stream.  These  would  appear  as  mathematical  analogous  forms  to  serum 
thrombotic  accelerator  activity,  which  appears  to  enhance  thrombus 
formation  in  the  flowing  blood.  An  attempt  could  then  be  made  con- 
currently to  change  channel  boundary  conditions  to  conform  to  the 
periodic  pulsatile  mode  with  viscoelastic  conditions  present. 

Ancillary  to  the  theoretical  studies  has  been  the  initiation  of 
experiments  for  in  vitro  verification  of  presently  evolving  theoretical 
flow  studies.  A mock,  circulation  apparatus,  bolstered  by  the  addition 
of  photostress  equipment,  is  presently  being  employed  to  evaluate  fluid 
flow  motion  in  the  transval vular  areas  of  the  heart.  With  the  inclusion 
of  another  section,  it  will  also  be  used  to  study  atheromatosis  sites 
in  vascular  channels. 

There  are  certainly  weaknesses  present  in  the  system.  Superposition 
of  a curvilinear  boundary  onto  a rectangular  finite  difference  net  is  ^ 
tenuous.  The  truncation  error  may  be  enlarged  and  the  presence  of  the 
non-symmetrical  finite  difference  stencil  can  lead  to  additional  error 
si  well  as  possible  propagation  of  the  round-off  error.  Many  flow 
configurations  are  complicated,  such  as  the  combined  renal  and  abdominal 
aortic  bifurcation  sections,  which  result  in  many  mesh  nodes  typed  as 
"exterior"  to  the  solution  space.  Computer  storage  must  then  be 
allocated  for  these  "unused"  nodes.  Consequently,  sufficient  detail 
■ .innot  be  allocated  to  areas  of  flow  where  detail  is  required  and 
such  problems  as  "flexing"  boundaries  compound  the  difficulties. 

uch  weaknesses  are  being  attacked  by  now  computer  "hardware"  studies, 
in  the  case  of  storage  problems  and  "software"  studies  in  all  situations. 
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The  storage  is  presently  being  tripled  in  capacity  which  should  allow  a 
finer  mesh  for  greater  detailed  study  at.  a specific  point  without  losing 
gross  evaluation  at  point  distal  to  the  obstacle.  In  turn,  another  view 
of  the  problem  of  a curved  boundary  has  been  investigated  and  it  remains 
to  be  verified. 

Another  weakness  of  the  presently  discussed  technique  is  that  it 
employs  vorticity  and  stream  function  values  as  the  primary  dependent 
variables  within  the  computational  mesh.  The  advantages  were  that  only 
two  dependent  variables  were  needed  and  mass  was  easily  conserved,  since 
the  use  of  a stream  function  identically  satisfied  the  mass  conservation 
condition.  However,  the  disadvantages  associated  with  the  vorticity  - 
stream  function  method  are  its  poor  extension  to  three  dimensions  which 
requires  the  introduction  of  a vector  potential  for  velocity,  and  the 
difficulty  encountered  in  satisfying  free  surface  boundary  conditions. 
Also,  true  channel  flow  remains  to  be  described  and  pulsating  motion  and 
viscoelastic  concepts  are  difficult  to  investigate  with  the  present  nodal 
parameters.  To  overcome  the  difficulties,  the  Marker  and  Cell  method  ( 
is  being  intensely  pursued  by  the  second  author  of  this  report.  This 
work  will  be  the  bases  of  the  second  report  in  this  series.  The  method 
uses  pressure  and  velocity  as  primary  dependent  variables  and  allows  free 
surfaces;  the  latter  useful  for  pulsing  and  viscoelastic  effects.  It 
also  maintains  accuracy  with  a minimum  of  computer  time  through  the  use 
of  a corrective  procedure. 

The  computer  experiment  is  basically  designed  to  supplement  other 
techniques  and  to  allow  new  fluid  flow  investigations.  The  contribution 
of  the  computer  graphics  approach  is  the  ability  to  present  solutions  to 
certain  previously  intractable  situations  and  in  detail.  Hopefully,  in 
the  future,  solutions  to  non-steady  flow  problems  will  be  in  the  forms  of 
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a computer  directed  optical  display  continuum,  copied  by  movie  camera,  with 
the  allowance  of  desired  changes  in  flow  configuration  upon  call  from 
computer  programming  storage  banks. 

In  the  future,  usage  will  be  made  of  the  head-mounted  display, 

64 

developed  by  Sutherland  , and  presently  being  augmented  at  the  University 
of  Utah.  Two  miniature  cathode-ray  tubes  are  mounted  on  the  user's  head; 
one  tube  on  each  side  of  an  eye  so  that  prisms  then  direct  and  form  a 
stereoscopic  view  in  front  of  the  user's  viewpoint.  As  the  user  turns 
his  head,  the  theoretical  perspective  object  appears  to  remain  stationary 
as  the  viewer  moves  into  an  appropriate  viewing  position  to  examine  an 
object's  particular  feature.  It  therefore  appears  feasible  for  one  to 
project  himself  into  a computer  displayed  prosthetic  heart  valve  and 
watch  "theoretical"  blood  rush  towards  him  from  the  atrium  or  turn 
around  and  view  the  valve  occluder  moving  away  or  towards  him.  Feasible 
too,  would  be  viewing  of  flow  about  an  atherosclerotic  lesion. 
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■n  FOR  XYMAIN, XYMAIN 

DIMENSION  SN (61,31) , SO(61,31),  VN(61,31),  VO(61,31), 

1 TA (61,31) , I FP ( 2 ) , 

2 BND (61,31) , A (61 , 31) , TB(61,31) 

DIMENSION  SI (61, 31) , S2(61,31),  S3(61,31),  S4(61,31) 

DIMENSION  BPF(100,5),  BPS(100,5),  IBF  ( 100 , 5)  , JBF  (100,5)  , 

1 IBS ( 100 , 5) , JBS (100,5) , NBPF (6) , NBPS(6) 

DIMENSION  l VS  (250)  , IVB(250),  .IVB(250) 

DIMENSION  ISB ( 250) , JSB(250),  JSE(250) 

DIMENSION  ISB1 (250) , JSB1(250),  ISE1(250) 

DIMENSION  XBP ( 100, 5) , YBP(100,5),  BLN ( 5) 

DIMENSION  EV (60) , ES (60) , ITERS(60) 

DIMENSION  BP1 (5) 

COMMON /COM 3/  ISB,  JSB,  JSE,  KSB , ISBl , ISEl,  JSB1,  KSB1 
COMMON /COM 1 / BND,  M,  N,  Ml,  Nl,  M2,  N2,  M3,  N3 
COMMON/COM2/  A 

INTEGER  BND,  BLANK,  SWITCH,  VTOBIG,  DISPLY,  EDIS,  PN , PNCH , 

1 STEADY,  S STATE , REITER 
LOGICAL  RESTRT , REDO,  EPLOT , 

COMMON/INTRC/  OMEGA,  R,  NFRAME , SFVBND,  SFVOP.T,  BPR1 , BPR2 , 

1 BPR3 , BPR4 , DISPLY,  DR,  DT,  EVOR,  ESTR , IWRITE, 

2 NIPR , RESTRT,  REDO,  EPLOT,  SWAPON 

COMMON  RELAY/  IFP,  NFILE , XBP , YBP , BLN , BPF , BPS , I BF , JBF , I BS , TBS ,NBPF , NBPS 

DATA  BLANK/6H  / 

DATA  PN/'PN'/ 

DATA  STEADY/' SS ’ / 

I TTIME  = EXTIME (1) 

CALL  DATA (M,N,END,A) 


READ  1,  R,  DR,  DT,  NFRAME,  IWRITE,  IPRS,  DISPLY,  SSTATE , PNCH, 

1 BPR1 , BPR2,  BPR3 , BPR4 , OMEGA,  DOMEGA, 

2 SFVBND,  SFVBND1 , SFVORT 

1 FORMAT (3F10. 5,  415,  4X,  4A2,  I2/4(5X,  F5.1),  6X,  F5.2,  7X,  F5.3/ 

1 3 ( 7X , F5 . 2) ) 

IF (DISPLY. EQ.O)  GO  TO  2 

CALL  RELOAD 

CALL  TYPLIN( 'READY tiA' ) 

2 NIPR  = 1 
CALL  NOSKIP 
EPLOT  = .TRUE. 

PSCALE  = 1.0 
J.7MAX  = 60 
MSP  = 1 

NSP  = 1 

I F ( N . GT .31)  NSP  = 2 

Nl  = N-l 

N2  = N-2 

N 3 = N-3 

Ml  = M- I 

M2  = M - 2 

M3  = M - 3 

H = 1 . /FLOAT (Nl ) 


H2  = H*H 

DH  = FLOAT (Nl) 

EVOR  = . 0750*DR 
ESTR  = 0 . 000050*DR 
DTP  = 0.0001 
DTI  = 2 . *H2/DT 
DTDX  = DT/H2 
IPSKIP  = 0 

DO  3 1=2, Ml 
DO  3 J-2 , N 1 

S1(I,J)  = 1 ./ (A(I+1,J) * (A ( 1+1 , J)  + A (1-1 , J) ) ) 

52  ( I , J)  = 1 . / (A(I-1,J)*(A(I  + 1,J)  + A ( I- 1 , J) ) ) 

53  (I , J)  = 1 . / (A  ( I ,J+1) * (A(I,J+1)  + A ( I , J-l ) ) ) 

3 S4(I,J)  = l./(A(I,J-l)*(A(I,J+l)  + A (I , J-l) ) ) 

4 DO  5 J=1,N 
DO  5 1=1, M 

SN ( I , J)  = 0.0 
SO(I,J)  = 0.0 
VN ( I , J ) = 0.0 
^A(I,J)  = 0.0 
TB (I , J)  = 0.0 

5 VO(I,J)  = 0.0 
RESTRT  = .FALSE. 

REDO  = .FALSE. 

BP1(1)  = BPR1 
BP1  (2)  = BPR2 
BP1  (3)  = BPR3 
BP1 (4)  = BPR4 
BP1  ( 5)  = BPR5 

C STEADY  STATE  DATA 

I F(SSTATE.NE. STEADY)  GO  TO  9 

READ  6,  ( (SN (I , J) , VN (I , J) , I = 1M)  , J=1,N) 

6 FORMAT ( 9F8. 1) 

7 DO  8 1=1, M 
DO  8 J=1,N 
SO(I,J)  = SN  ( I , J) 

TA  ( I , J ) = SN ( I , J) 

VO(I,J)  = VN (I , J) 

8 TB  (I , J)  = VN  (I , J) 

C UNIT  14  IS  A FAST  DRUM... USED  FOR  TEMPORARY  STORAGE 

9 REWIND  14 

WRITE (14)  TA,  TB, 


CALL  VBND1 ( I VS , IVB , JVB , KVB) 
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Mvj.i  t)»  ■.,#  ;n  » t ax#  srvoPT#  sfvb;io»  kvor.  r.tr# 

1 ...  •>.#  Ml  hit  , T|,  , jF  VI  i:  1*  OMEGA#  I'Of  Kr  A 

)1  (•  o v -Li  ( 1H;  , //dux.  iOnf  o ,RAM  PAt,  AMf/TLi 'S//20X  » 2HN=I.V/20X»  2MM=I.S 
\ / / < > a # ,’iv;f  1 Ul  =l-  . G//20  X . GHT’MA  X=F9 , X » 

•St-y.Kr-*  ( 1 ■♦.,*//<:  C'a  r ♦ ,1  VO  .i  F14.7//2.)Y»  GHFVOW-  F14.7//20X# 

' b-i..'  T‘  14  .?//’)>  * UJPAX::  I4//20X.  2HH=  Fll.7//  POX# 

4 4l  |-=,  114.///  ,0a,  • I/hL-’.  F 1 4 • 7 / / 2<)X»  * SFHP 1 = • # F14.7# 

5 //.  X.  'uPl^A-'#  U'..?#  //POX.  *OGVLGA='#  F14.7#  1H1> 

IF  l iu  t G o , f • I ! # (j  L 1 1 j)  2 . rc  It 
hi...  1./  4 

4 l mi  ■ ( 4 ) •'!  • I . » I «F  IlL  * T # T M A a # t<  • OLL/I  iR 

OWGHi-  - J*  I 

(..ALL  IN'  JT  ( J , . J,N  v,:  !,.) 

CALL  1 NOUl  11#  + # M # i'.  V / \ ) 

F L mu  l 4 ) i i )F  » .41  F # .l),-’1 

L'O  i c I - 1 # i iF  1 LL 

h L At  v ) iF  » ( x !jp  ( K # I 1 # Y , F ' ( K # I ) # K — 1 » IK) 

U LL'41)  r i,; 

SMlL  = i 

ip:,«  = 1 

I.T  I KL  - 0 

hi  - Ltd 

I ^ u . u 
IF  - T . !) 

l 0 12)  a - 1 # i .Gr  I P # /i 

lF’Shil  r Il'GKlP  4 j 

If  ( 1PGF  IP.V.L,  if  i < S ) IlV-'.Tr’  = 'j 

I T l t.  - hi  I f>'i  ♦ 

1-14  i.0*Lil 

IF  ~ IF  4 2.U‘(JL  V-i'/h 

fi  - • i ♦ .5 . l *uu  jl'f; 

C A l l J t/’U  1 ( .1  # 4 # 7 A # Gir.  ui  ) 

C All  1 Nf  U I ( 1 » 4 » I i » i)'<  ol  ) 
l-'F  „i  / 14 

a h 1 1 1_  ( 1 4 ) 1 A.  1 it  I l i if. 

C all  1 1 JOU  1 ( 1 # 4 # Yj  # .agi  ) 

call  mcu i u # 4, v'- ) # dvoi  ) 

LAll  I Fir  U I ( 1 #•*  #S4#.iwj!  ) 

L ALi-  JilGUT  l 1 #4,  J hhv.Jl  ) 

1 J lOM  1 . .1  IF 

IF  (f  c r . ' A. . f:  I ) t r K 1 
t All  i u7 f lu 

Fl-.lUl  14 

14  f Gih-  >,  ' < 1 Mu  » tLX#  * ....  1 ill  1 AL  GUFbb *#///) 

C A . F iwr  (H#  T#f  'AX,<',L1.JF>,NGM»Mi»VU,L#  I T’  H » JJ.  T # IP  #h»;0#M#N» 

1 f-  1 # . 1 » I f # 1 f # JhF»  i.'  r’F  ) 

IF  ( G j7  A 7 1 . ; ,F  . S 1 1.  A l>  7 ) « . 7 0 <r 0 

i jO  l „ 1 r 1 » 

LG  * o J — 1 # !! 

GO  ( I » u ) = of  i ( I » J ) 

7 A ( I , J)  = S!l(  I # J) 

V U ( 1 » J ) = v 1 1 ( i » J 1 

lb  It  l 1 # J)  - VI  l 1 #J) 

hF  ..  I . *1  14 

t (•  1 | ( ( 1 4 ) T A » T : t 
GG  lo  10 


! C > M 


17 
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CO  to  ( 1 7.  il  ) , 1-.  ill. 

K fc. : i , * it 

K I 1 L-  (4)  M»  f . iFILcr  [ I • I 1AX  * l'i  OP 
UmOi  .u  — i >1  # ’■ . 

call  INCUT  ( 0 » » I ..,.N„r  r) 
call  iMOUl  (U»4f  fcii-.v.,,;!  ) 
itK  i I c.  ( 4 ) 1 [if  t Ji  ll-  * : . 0 > 

LC  Id  I ~ 1 • it  iLl 
IK  = i U'U) 

1 0 while.  ((,)  I y , (a  i (K,’)»  yl;F  ( k » I ) » K - 1 » l K ) 

lc  bill  r x 

lh’L.  ■ - 1 
1 ~ o . 0 

It  = .0 

i.f  i;  _ - j 

It  ( d-jT  >.  Tt  . M , ‘ .7 1 . ,:i)  Id  l,j 

r;Ti,,u  - 1 

It  — it 

1 _ r<  * QTP/i H 
L T 1 - 
l I - I TP 

IN  i l .0 
CO  lu  07 

*_♦  **♦♦.***  + **,  , 4 * * » ,**»****4i*  ******** 

o o oo  — j 

Kt  i TU'  - O 

if  ( j i ; Af'r.Lc.s^Tn  r i • - 1 tlp  - j 

n\.D'i  - U.II 
I I i/i/i .i  ~ Ij  • 0 

nvoKi  - u.o 
T T bT t\  s U • u 
I T • 1 jC  = U.O 
T - I 4 OT 
If  ; IF  ♦ uI'/k 
f.Tl  -i.  x 1 1 i t t + 1 
t " I - 0 . b 

It  (MIMt  .Li  .l,H<«' ) uO  Ij  07 

v O T ^ ( ?1 ic4 I > > * I T 0 1 1 

«^1  bl  * - 2. *t ,<‘/d I P 
bf  = Tn 

t _ i ♦ o.*i  rij 

f .t  I-  A I*'  l ~ s.ii 
L-A  - f ,n 

It  = i . 

S* 1 I U’  - c 

It  (t-.j.  T‘  .1  , Cil)  dO  Id  . 

2t  call  ckp!:.  j 

PUi.Cti  u.  ( ( 111.  ( 1 » J ) » 4\  (1»J)»  I = 1»M)»  0=1  »N) 

CA,.l  CKPIwt 
CALL  ‘ XT  r 


1< 


24  GO  TO  (25,26),  IWRITE 

26  CALL  NT RAN (4,9, 11) 

2o  CALI-  EXIT 

27  DO  29  KK=1 , NFILE 
K2  = NBPS(KK) 

IF(K2,EQ.O)  GO  TO  29 
DO  28  K=1,K2 

I = IBS  (K,  KK) 

J = J BS  ( K , KK ) 

28  SN ( I , J ) = BPS (K, KK) *TF*BP1 (KK) *R 

29  CONTINUE 

IF (NT1ME. EQ. 1)  CALLOUTFLO 
DO  31  KK= 1 , NFILE 
K.1  = NBPF(KK) 

IF(Kl.EQ.O) , GO  TO  31 
DO  30  K=  1 , K1 
I = IBF (K, KK) 

J = JBF (K, KK) 

30  SN  ( I , J)  = BPF (K, KK) *TF*R 

31  CONTINUE 

IF (REITER. GT.O)  GO  TO  37 

^****************************************** 

C COMPARISON  OF  TIME  ADVANCE  BY  ADAMS-BASHFORTH  4-POINT  FORMULA 

C AND  QUADRATIC  EXTRAPOLATION  

****************************************** 


32  REWIND  14 

READ (14)  TA,  TB 

DO  33  1=2, Ml 
DO  33  J=2, Nl 

IF (BND ( I , J) . GE. 1)  GO  TO  33 
TC  = TA ( I , J) 

TD  = TB ( I , J ) 

TA ( I , J ) = SO(I,J) 

TB  ( I , J ) = VO(I,J) 

SO(I,J)  = SN (1,0) 

VO(I,J)  = VN ( I , J) 

GO  TO  (34,33) , IPSW 

34  SN ( I , J ) = 3.0*S0(I,J)  - 3 . 0*TA ( I , J)  + TC 
VN ( I , J)  = 3 . 0*VO ( I , J)  - 3 . 0*TB ( I , J)  + TD 

33  CONTINUE 

37  REWIND  14 

WRITE (14)  TA,  TB 
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1 1 ' lv*li)  T t_ ' 1 . • • 


\CC  PY  Ai  A viS*BASt 'F ORTH  OK.Y 


Tm  7Y, 


) L 1 1 


lflji-*  ( I > o*l  ) *l  • i *wv«  j * Y .OR. 

’ ,)IY  ( I *1  * J)  .L.  ..i  .>J-.  •':!'(  1-1  * .j  ) , t_  o . 1 ) V>  10  3S 


FN'i  - ( '•  0 < x , O ♦ 1 ) 4 V .,  ( t » J-l  ) -t  '•/(■’  ( I ♦ 1 » J)  + VO(I-l»J)  - 4.*VC(I»J) 

1 - ii.t" * i ^ lit.  i j)  - voci-i  >j))*(so(i'j+n  - sr  (i.j-n) 

2 4 I,.,.  )♦(,/,(  I ,J4  1 ) - 1/tMliJ-l)  ) * « SO  ( I ♦ 1 »J)  - S0( 1-1 »J> > )/H? 

4<U  _ (ll)(ilu4l)  4 1,.  (1,0-1)  4 It  (I  + l.J)  4 Tll(I-ltJ)  - 4.*TMI#J) 

1 - >.<  )*i  I <1  ( 1 4 I .o)  - TLtI-1  tJ)  )*(TMIfJ+l)  - TA(I*J-1)> 

? ♦ (l,., 4j)  - run  »j-d  )*(TAt  i4i»j)  - i At  i-i  *j) ) )/h? 


tb*.  = (TJ(I.j4i)  4 I (..J-I)  4 I (I41,J)  4 T J ( I ~ t f J ) " 4 . ♦ TT  ( I » J ) 

- 0 • <JS*  ( I..  (14,  fvj)  - Tot  I-l»  J)  )*  (1C  ( I * J+l ) - TCtl.J-l)) 

1 4 j./s+t  f^(l».i4]  ) - fL(  I» J-l)  )*(TC(I+l»J)  - TCtI-1  »j)  > )/M2 


v f • ( 1 # _ ) = \jC  ( I * J ) 4 (,  . ,i  *1  u j - 16.1)*FM  4 S.0*F*J2)*l)T/12.0 

Pv  TO  3r 

3b  VN(1»J>  = 2>..,»V0(  1 ...  ) - .'••(.)*  Tl  (IrJ)  + TD  ( I » J ) 

3F  CO  .Til  UF 

It.  .,iu  13 
»<  K 1 T l l 1 ^ ) v/l . 

M\  PR  _ *TF*PSCALk 

IF  1 <F  1 Lt  • Cl  • 3 ) 0 LL  _/i  1 lL 

iF  ( , i joT . bwAPOh)  PP  I jT  Nil!..!,  Kp 

-is  PCV'1/.  r ( 1Ho//3iJA*  P.'IC  -.CM  TtIWS  FOP  NT  I ML-  * Ir  » ?X»  2HR=F6.0/) 

Jl  ( . 1 4 i . O/i  AFOi  t ) PR  i . . T Hu 

44  4!  Fuu  • ),T  ( TUl  r i».  bn  IT  t;  , lUAr  1 JHVOP  T I C I T Y » 15X.  6HSTRkAN». 

1 4*1  ttl  jU'  lift,  .-I!., OK  PAf.  A 4,  * 7X  . ion  F <°0P  SUM»  3X» 

? • T bTU • » '4  X . *i  'J'.:  T * » 4 a » * TVL.iL1//) 

*♦♦**»♦♦*  » ♦ ♦»  4 .**4***»************** 


....  A I u LOOP  .... 

in  z i*]F  4 . b 

wlj  ~ ( 

t L 1 b - 1 
f.kL  - 0 

41  IF  (ji  .PL  Y , ! 0 . 1 ) .ALL  j 
11  1Kl'  T T ) G<;  to  4 
Ik  (Ki_l/0)  CuLL  Aunl  t 
10  Ht  1 Z 1 4 I* 

10  it  Jr  1 • i , 

T h ( I » J ) = blit 1 fj) 

4c  1 A ( I.J)  - VfJtl  »J) 


T <UP 


OJ  - JO  4 i 

T M I bO  - fc-  >.  1 I • 4 ( » ) 

C Al.L  ST i ti.P.t  V , T TF  n.LS.OI  1.  (>A  » L’Jl.w » JJ  . L jl'MT , A » Si  » SP » S3 1 S*4  • 
1 ^Ok'i_  jh  * O.  ’i  uAl»LSTK) 

TSTr,  < XT  IP'L  ( t ) 

CALL  ''u’.UtbU*  VI. » IV,  , lVb»KVl  ,'jFVIiNI  » h 2 » A , SFVHH1 » SI  * S2 . S3 » S4  ) 
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TV;  , : tXTl'iu(i) 

C AL  _ V lie  T ( Vf  , \i  0 * bi  i * , 

i v jk  i _ t/.i  r ,t_  ( i ) 

ITf.ljC  r TTf'i'jC  ♦ rii'i. 

I TVu’A  ~ T l Vl  If.  + 1 V i > i . ^ 

II  ,oK7  r TTVGkT  4 rVvjHT 
T T T t\  — T f j T k 4 1 j T K 

I T iMt.  — ( r ' ISC.  4 TVnMj  ♦ 
TTcOr.f  : irVoNU  ♦ Tlrf'1. 


H.SF  VORT  *!.IPt  ) 


rVOKl  4-  TSTR  ) / f»Ct . 4 
T 4 TTSTl, 


TTIMF 


lLj  - O.u 
L L V - O.u 

I O 4 j 1 — 1 * .*1 

L U *4  w Jr  1 9 • w 

It*  - Af'A\l  ( c — V » At>S(  T ( I * J)-VM  I *J)  ) ) 

‘O  t_t.j  - A'  A/.  1 ( t_c  j * r«t»S  i T , ( 1 » J ) — Sfj  ( 1 * J ) ) ) 

ui  <.b  I ~ 1 » . > 

00  *»b  J-l»u 

II  l L I4i  ( T * J ) . Jl  . 1 ) uO  T 1 . 4*4 
UU.J)  = /.I  j ( T A ( 1 * J ) -V  r i ( 1 4 0 ) ) 

bl/  TO  4 0 

4*4  T A ( 1 * v ) - u • U 

nb  (.001  i.1  Of 

IF (ElV  .1  T.lVOH  . \ t ).  I S.LT.FSTr^)  GO  TO  <+o 
VKjlf  r l 

IF  (ttV.oT.lu.tlG)  y/TuolG  - el 


il  ( . ■ f T . bfl  APGl  i ) 

1 PkI’  T 4,-,  jj,  LXW  t S#  ITc.R*  O^EGAl,  E'JEa*  TSTF' , TVORT.  TVBND 
4t-  F Ou*  AT  ( aH  * b X * I b » i p j . n * 3X » F20 • 6 * 3x  » lb*  SX , F 1 0 . 6 * bX  * 

1 * 1 2.b,  3FG.2) 


LV(jo)  - ; L\i 

irLRb(jj)  = itlr 

t b ( Jj  ) = Lt  r> 

11  Ill/1S,c.(j  .0 ) GU  TO  *4 7 

11  UiSPL  Y.EQ.l  . Al  .1/ • tPLOT)  CALL  ERRD  IS  ( NT  I ML  * JJ  * I TERS  * E V * ES » IR  * 

1 1 VOF\  »ESTr  »M*N»bLl<*  A.oP  * Y )P»iiF  ILL  * TTIME  * TVBND  * TVORT  * TSTR  * T^ISC  * TA  * 

? TTLOMP.  TTM1SC*  TTVO.a  *TTV(jNl  * T1  STR  »NFL ) 

47  GO  l j ( 4 3 » 4 y ) » vTOuiij 
4 is  IF ( JJ , LT • JJV AX ) GU  To  ■*  1 
CAlL  r XI  T 

4y  IF  lb*APON)  ou  Tu  t>4 

I'RIiJl  4 0 * Jj*  tLV*  t.r  St  ITLK,  OMl  GA 1 * F ME  ft  * TSTR*  TVORT*  T VP  NO 
PR1.J  bO*  TT  jTR  * TTVr.f  , TTVOKT.  T1C0MP 
bu  FORMAT  (lltU*  2JX*  * TO  f AL  TIME  (SEC)’/  24X*  *STREAW=**  F7.2/10X, 

1 »VolTICirr  (BOUNDARY )r* f F7.2/10X*  ’VORTICITY  (INTFRIOR)=»  E7.2, 

2 /VoX*  * TOT AL=» * 17.2) 

1*  t .«  1 . ,1  13 

Ft  AO ( 3 3 ) ltt 
TLn  - il.U 
00  b 1 1=2  * >1 
Go  b 1 J ~ 2 * i J 1 

IF  (;ti«  ( T * j)  .GL.O)  jO  I / jl 

IF  Ut  .Gl  . vH'jf  VUl  1 * J)-l  A ( 1 * J)  ) ) jO  To  bl 

tl-'i  - at  ■af  vr  * c i *-»)  - i/.<r*jn 

IV  L I 
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J V - V. 

bi  c c . , r i u 

► Kil.l  )"»  1 A ( I V » JV  ) » ' (1  V,JV),  IF  U IV*  JV 

1 Oh,.,Al  ( 1 Hu  * ‘-a*  • C ILIi  (A. .-.)-**  F10.*/«!0X, 

1 *AbTUAL-*,  F 1 0 . _/ 1 . i > * * ' m X FlPf=».  * AT  POINT  ** 

r w . • , • * I o ) 

S3  1 1-  ( V f 0 i.i  I G .OT  . 1 ) C ALL  ' ] 1 

f PINT  ^ 1 lo  SIOKl  . ) A T M . . . . , 

b4  . 

IF  (i.riMf  ,LL.  1 ..*!•>.  h(  I TlK . fc  ) • G ) 30  T<.  So 
i P _iK  i > - 1 i’SK  1 1 ■ *■  1 
U ( lPbKlP.LT.IPh  .,)  vu  ! 0 5 b 
1 Pi<  I • * T r i 
A f"  jh  Ju  ’ ^ L 

bb  CALu  t X 1 1 1 1 ( l 1 . ■ 1 * I : X , 1 • , ibP  * m jP  , S'  I * V Fj  * I VvK  I TF  * T Tp  R , Jj , T * IP  * RN'D  *V  * M * 
1 Vi, M.  I PF  1 .'!  , [P,  I : , F * J - ; F * NBFF  * TO) 


CALU'LA  Tl  L X 1 1 Flo*  i OF  ILL... 

■jo  l all  < utFlo 

i.  i\  I T F VALULb  OF  bl  I A-  AN,  V OPTIC  I TY  FUNCTION  OM  JM 

TO  SAVl  IN  LmSI  Kl.“  initialization  is  jesikfd.  . . , 

K E * i < «[  1b 

»■  P I T L ( ll>)  bIJ*  VI. 

( o oO  - . F,i  LSL  . 

Ft  1 1 l = i\i  Iltl  + j 

1 F l * . I I ,*lfc  . L l . i . 1 • R .Lt . 1 ) GO  To  ? 7 

V>0  TO  rJ'l 


♦ ».♦,  *♦*♦»♦♦*  + *♦****♦***** 

bUbr  01  T 1 Jl  AO  A i i • 

TFiIS  HbjTlJt  l-Ii.IT]  .LI^L'j  (if  bT‘  t’Ai^  VJJ  VORTICITY 

MJI  Cl  luNS  10  T . i(  Ik  . A lUF  b AT  THE  MLGINNINj  OF  T HP  TIME  STEP 

FT  a I ni  1b 

RL  Al  lib)  b :»  vh 

JO  — f. 

KEl'O  ~ . FMl.  Si_» 

»•'  t T OKI 


bUriKOl  T I IJL  Pc.RTi<i 

PLk  l = Mil  I 

IFlKI  = If  l F I 4 1 

IF  l 1PI  NT  .oT  .1 ) nr  TUHl. 

VNi ma  s O.o 
1-0  1 1 = 1 n- 
LiO  1 J=1 

IF  (Fiji  ( f , j)  ,;,E,u)  oO  T ' j 
VlMAA  = A,.  AX  1 ( Vu-'AX  , (I,J)) 

1 CO.il  i'  tjf 
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V I ( l r til)  - 1 • S*  V ji1  A X 

VN  tl  / » 11  > - 1 .b*Vi  V.Aa 
Fit  f uK1 


v. 


c 


LMlhUl  1 I i'll.  0UTFLO 

• •••I  OviLK  « • a • 

If  ( !.F  IL f .l?  .0)  GO  1C 
1 = . 4i  PF-  ( J ) 
i F ( i . I (,.U)  GO  10  4 
ji  = inr(i.i) 

12  = IDF ( I ,3) 

if  uii  imf  .oT.i  .ok.  uu.  T.i)  00  ro  2 

Ft,  1 - SlMd  1 ,2) 

TL ' r - 4(  i 1 ,3) 

U - (Slid  2*2)  - GM  ’ I * 2 ) ) / FLOAT  ( 12-11  ) 

LG  1 J=U»I2 

S0d»J)  = TLM2  F TLM*'  i.OATl  1-11 ) 

1 SMI*;)  = TLM1  + Tt  M*F  LOAT  ( I-1 1 ) 

2 K = u 

t 0 3 ! = 1 1 , 1 2 

K - K f 1 
bull*!)  - bM  1 ,2) 
i IT  (K.  3)  = SM  i , l ) / dF  *Fi ) 

J VM  If  1 ) = V N C I ,2) 

. . . . I K . . . . 

4 IF  (fjMU  • L T . 4 ) GO  TO  r 

I - hi-PF  (4) 

IF  ( i .L  L.U  ) 00  To  u 

II  - 1 1 IF  ( a » 4 ) 

\C  - IHHIrUl 

If  ( i . 1 I J'T  . GT  . 1 .OK.  Jo.., T.I)  oO  TO  6 
Tl  N 1 r SM  1 1 * Ml ) 

U.i-12  = SN(  11.U2) 

I tf  = (SMI2»01)  - SM  11  ,N1)  )/FL0AT(I2-ll) 
00  b I = U,12 

bin  1,1  2)  =■  TLM2  *■  TL.'l*FLOAT(I-U) 
b SMI,rl)  = T + Tt,.imOATII-Il) 


U l*>  = 0 

DO  7 1=11,12 
K = K + 1 

Slid,!,)  = GNU  nil) 

L,F‘F  l tv,  4)  = SMl,N)/ll)  «f< ) 

7 VMl»i  ) = VMI»Nl) 

....i  I'jbT*,.. 
h 1 = ui  pr  (2) 
vj  1 — Jilf  (li2) 

J2  = JlT  (1,2) 

IF  ( ( « I I -If  . oT « 1 .07.  Jo.uT.  1 ) oO  10  IF) 

1 EDI  = SN(,-1,J1) 

TL,-,  = <SN(M»J2)  - GDC  ,J1)  ) /FLOAT  (J2-J1) 

10  y J= j 1 , u2 

SMM1»J)  = T F M 1 ♦ TL,-i*rLuAT  ( J-Ul  ) 
y SU(M2,J)  = Tt.Ml  1l.-,*FLjAT  ( J-Jl  ) 

1U  WH  = f - 2 
00  11  U=J1,J2 
DO  11  I , Ml 

sm  i ♦ i , j)  = sin  i , j) 

ii  v/iiin  , j>  = v n i , J) 

K = U 
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LO  Jr  J l # J J 

K - r\  f 1 

1 £ PPp  ( i\  * d ) - Si . ( * J ) / ( i »»'  ) 

• • • L Li  i . • . 

IF  (Ptl;T(  ) la l.L  PLF  I K, . 

Afc  Tijlv 

♦ ♦♦♦♦*»♦♦♦.,«.*♦***♦>*♦***♦******** 

bUuHol 1 T 1 N|  O l! b T L . 

M - « 

1 K1  =.  ( 1 ♦ 1 

if  (rXl.Gl.tji  1 i i ) P T Ok ’ 

Cl"  O I (.’i\  1 [■.  I r .<  1 0i\  n;  ,rACl  i L. ........ 

L - uL  '1  ( K l ) 

IF  ( li  I x ( A,  P ( l ) ) . L-j . 1 .OK.  I F I X ( XuF'  ( L t K 1 ) ) . t Q • M .OR. 

1 i »-  1 < ( 1 iM  1 ik  1 ) ) >i  « . . . )R  . I F 1 a ( YUP  ( L t K 1 ) ) . rQ  • ' ! ) Gu  T 0 ] 

K r.  1 
lKt.F  - 10uU 
JJ  = i Li  (r. ) 

l>C  F » -1  iJJ 

d I K t_  F z 1 1 . ij  ( i K c F » lii.  (X  >P(J.K)  ) ) 

TL-.  - O.U 
JKLF  - 0 
JJi  — JJ  — 1 
L’O  J L— 1 I w J 1 

If  ( 1 r 1 X ( Ar  P ( J » \ ) ) . , H.  . . r r ) O J K J 

I L - T 1 A *•  \ . 

Jht  I-  - JKl  F ♦ li  IaIi  iPlwiKH 
5 (-0  il  i.’  Uf 

1 t.,-1 1 - F Lu>\  T ( JKlF  ) / Tl 
jRt.Fi  = rt,-,i 

1F(Tl'  1-FloAI  ( JPLPl) .vl  . J)  GJ  10  4 
JRLF  t ~ Ji'tF  1 
CO  ru  b 

4 «JRLF  ~ Jim.F1  ♦ l 

b if  tnr  i mi  .gt  . mo ) '.->0  ic  t 
bilKLF  1 - i»i4  ( 1 » JRLF*  1 ) 

StiRLF  <:  - o il  1 1 JPt  F c' ) 

(0  fu  u 
O U I r\  t_  F z 1.0 
00  7 J-|  i jJ 

7 L.l,<tF  - A.-i  X r 1 1 (OI’<„F  i , b(XHP(J»K)-FLOAT(lRrF))) 

1F(^1-f;F  .Ll  .U.OOJ1  > aP  F = lMtF-l 
b = A(  II-’tF  ♦ 1 »JKF.F  1 ) 

jNk.F  ] ^ .b*ii<(i.+  ) *r.  ■ ( 1RLF  -2 # JKEF  1 ) - 0 * ( ? . ) *SN  ( IRE  F - 1 * JRFF  1 ) 

1 ♦ . u*  (!.♦  i>  )♦(,.♦,)*  J l I <f  F » JREF  1 ) 

L>  — A ( 1 K e.F  ♦ 1 » JF  t_  F .•  ) 

bl  iKi_F  <•  r .•*■»*(  1 . ♦ i»  { 1 1 ’F.F  -<•'  # JRtF  c ) -ll*{..+U)  *C>N  ( I PE.F- 1 » JRF  F?  ) 
1 ♦ .b*  ( 1 . ♦ .»)  ♦ (2.  ♦ •)  * j ( I ,<FF » jRFF  ) 
d II  z lOliU 
U r n 
jl  = iono 

d*  - Cl 
CO  ')  ,:|»oJ 

II  = If  U(  . 1 . IF  lx(  A, I ( J » K ) ) ) 
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12 

r 

AV0( 12* 

ir  1 X ( X ,! 

( J 

* K ) ) ) 

Jl 

r 

InG  ( Jl  * 

IF  IX  ( Y.J 

( J 

* K ) ) ) 

9 J2 

— 

A X 0 ( J2 * 

IF  IaIY.s! 

( J 

* K ) ) ) 

il 

— 

:i  - i 

12 

~ 

12  ♦ 1 

Jl 

— 

Jl  - 1 

J2 

z 

o2  ♦ l 

l U 

iU 

1=11*  12 

uO 

10 

J”  Jl  * J tl 

IF 

t L#»  i 

' U»v).NIl 

• 1 ) Jt«  1 

J 

10 

bU 

(1# 

J)  = .bMSIWEFl  + 

S 

111  to 

‘Ji 

.JF 

GO 

To 

1 

v 


i 


j 
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I F uK*  t •'  1 1 1 , t K 1 i T 


00'  i<  j T i 'll 

f •'  1 .M  1 (O.W 

I t T 1AX  f\<9  * 

,r  > , JbP 

9 ‘if  4 9 J 9 

l 4 R J 1 ( , lb, JJ,T, IP, 

f .V  , 

P I'i  9 M 1 9 . 1 9 

1 PR  I r.  T , TF  , I ‘ >F 

* JPF  ,f.'(PF  , TA) 

I l'-u 

i 1 ' f (1oJ,M)» 

Jhr  ( 1 Uu  p a 

) , NoOi 

{<  ) 9 T 

A (61  * 26 ) 

1 L'O  I Oi 4 
I!  Tt  Ol  S . > 

1 . ( 1 1 ,/.*,)  * 

ij  9 iLA  W. 

( • > 1 9 C t j ) 9 

■ •NL'(ol,?b)  , 

K.OUT  ( f 1 ) 

i A Tp.  i L ’'it  /t>i  / 

MdilKj  ~ ^*,4 

00  1 u ( 1 * t) ) » If  . < i 1 1 
IP-  i(.l\Lc(SNi'l'Mi) 

P P 1 f 4 1 2 , 1 , *J  U , r ' , I , K 

L f Of.  ,mT  UMUf  *?*•  tit  U.j  \ JE  AM  F (JLC  f 1 0 >i  V .LUES//-  7Xt  JSHAV.  ST KF Af*  ITFN. 
1 Ijr  MM  vOfl.  :TrH_];,  3h  f - T ' • . 2 » M-  T=FO . b/ 57  X * 2Hw=Ffl.O) 

1 l h I = 1 9 • 9 MSM 

L'O  b jri  9 1 9 M jP 
r 1 - i *C,N  ( 1 9 J ) *■  . i 

a Oh  1 i J ) - i NmCD  ( i\  1 ) 

b If  I.H  ( I » o ) . t.w  . 2 ) KOjT(J)  - I ■ L r UK 
«4  PRIM  b,  IkOoKj).  u-l*N*NSP) 
b FOK\,mT(  Im  » bx#  blA-*) 

t>  00  T J ( 7 9 ;i  ) 9 1 U I f L 

7 '*NI1ii(4>  r 

0 All  I f Jf'0 1 ( o 9 + 9 iy  i * 1 1,. u < ) 

to  60  lu  («»M)  9 lr>Kl  IT 

y vh(i»  i ) - vM2*i) 

V N ( 1 * l ) — v f ■ ( l » 1 4 ) 
f = :■>(  Al  t (VM  ^ * t J ) 

f R I .N  r 10*  0 

10  F Cf . mHHIui  67a*  1 lo  T I C I I Y FUNCTION*  5X9  2HP-F9 .?.//) 

^0  lL  lrl,  9 >P 
DO  1 i J=1 , 4*.JSf 

K 1 - 1 * '/N  ( I * J ) -9  ,‘j 

KOUT  ( J)  = UoCt  (si  ) 

11  1 F ( tj I <4 [ ( T * sj  ) # c. ' j • t 1 ) is 0 ■ ) T ( %J ) — iLAf  >* 

It  PRIM  b * IKOUIM)*  J-l.N*NbP) 
lb  oO  Tu  (14,1b)*  I /Kilt 
14  CALL  INOUT  (D*4, 1/IJ,f4Wu-  ) 

lb  oO  I u (1002)  * 1 ‘ 1 1 . r 

l (>  00  1 o Ulc' 

U 0 ltj  J ~ 1 * ; 4 

If  (ij.t.  ( T * J)  .LOO)  i»0  f 17 

1 A (1  ) - II.L) 

oU  To  IP 

17  1A(1,J)  - ab/LA*  (SNl  1 * Jfl  ) - S,N(t*J-l>) 

1 b COM  i*  u* 

C A l L 1 Mf  ’ A i . 1 ( *0*  '1*  l,i  I-  * Jtj  f , i . . 1 1 F , 1 A ) 

F - j(.  A l L ( I A * W , , ) 
f Mil,!  1L  * V 

19  FOHMAT  ( 1H0»  J/X9  IvHu  Vf  LOCI  TV  FUNCTION,  5<t  ?HP=F9.2//) 

LO  LI  Ir  1 , ' , f<  SF' 

00  t u 1 , i i » * < .>( 

K 1 = • * T A ( I , J ) 

KOljI(a')  - ll4l)CU(-.l) 

?0  IF  ( (,  < I *o)  .1  >..i)  KOu  I l J)  = 


(LA:  K 


IK0.J1  ( Jl,  J^!  , Od.SP) 
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J • «-  1 


LO 

0 J 

I ~ l * 

1 V/ 

c~> 

.1=1  * i 

If 

( . . 

(I  * J ) . c 

/.,D 

1 A 

1 1 *v 

) = 0.0 

uC 

I J 

> 4 
C 

« ll 

Ta 

l i * V 

) = -U. 

'j  / . , X 

, J 

CO 

. w 

Ijf 

0 A 

1 )’  Ai  t ( M 

* N * 

f 

1 V 

Al  1.  ( i A * 

* i) 

PR 

i\l 

<0*  * P 

t >» 

f u 

* -1  /A  1 

( 1 H u * J 

/ >.  * 

i 'O 

C J 

1 = 1 *,• 

L 0 

J- 1 i . * 

■jp 

K 1 

• 1 

♦ 7 A ( l . v1 

) * 

MJ 

N 11, 

) - iltfj 

cou 

p s 

If 

( 14 

( fj  * w ) • C- 

4.c  ) 

e O 

f->. 

1.1 

'.i  t i r,00  I ( j ) 

DO 

cl_> 

1:1  , 

JO 

o 

Jri,„ 

if- 

( J.  A 

( l * ^ J • L 

lv  . (I  ) 

1 A 

(If. 

) — 0*o 

00 

lu 

? A 

i 7 

T 1 

.1  - 

.c/,  /♦  ( 

■ilill 

w. 

• 

. “'/..A*  ( 

1>N  ( l 

1 A 

( I * V- 

) = j0f<  t ( r i 

c.'U 

CO 

i ll  if 

01 

0 A L L i 

f 1 ' At\  T ( i 

» f 4 * 

l 

- ,( 

•H  c ( TA  ► 

* 4 ) 

f K 

i; . ( 

o.o*  P 

f Ui<  ' i\  I 

diii, » j 

7 X * 

00 

ji 

I r 1 * :•• . i'’ 

SO 

00 

Jo 

J- 1 * I-.  * 

op 

f 1 

— 1 

*T  A ( £ * w 

) 4 

r-  0 

U I I ,. 

) - INLj 

C 0 ( r\ 

JU 

if 

(ON, 

( I » o ) .0 

0.0  1 

Jl 

f 'i< 

INf 

.')  * ( K 0 0 

7 d 1 

JP 

f L 

r jh 

00 

JHU» 

T T f 10  I <-lu  A K Y 

1 i 

■ ■ C.  1 1 

ION  4|<r- 

<10  0 

uo 

J K 

K = 1 * 4 

L 

- 14 

PMkM 

u 

(L. 

0.0)  00 

TO 

00 

H K 

= 1 *L 

1 

- 1 

f-  ( X , 4K  ) 

J 

: vj 

f ( X * ISf1  ) 

If 

(o. 

> . 1 ) GO 

lu 

t A 

(1*1 

) - f A l i * «.  ) 

00 

T J 

<4 

1 

if- 

(c.' 

f .’4)  00 

Id 

Di-jl  H 4 


y 

■ -•  w * I 


I a < i » i » = r a ( i , , 1 1 ) 


GO  I j 4 

ihi.  t.D  w to  ^ 25 

TA(1»J)  = iA(<;»J) 

0U  i J 4 

J IF  ( i . L.M)  (,u  Tw  4 
7 A t #J)  = T fi  ( •'■"l » J) 

4 con  r ir  U( 

s a,i  if.jf 
Ht  T U l < ' 
t fJO 
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l l G‘  V .V.  • Vi  . ,L  1 

» U • ■ \y * 1 * j . V LilllijfiV'  tJi/btKV  ) 

i • * > J . V j 1 .*  0 ■ ) » i 1 ( .0  'J  ) > ov/  ( c.  0 (J  ) * i ;i)(ol  >?ti) 

i . 1 i * D ly  i*  J , u ) » wi  i_  i r (200) 

il  I i i ; i‘.  (.•  o>*  <2U0)  » ISFK200) 

lu1  j'  /Co.  .1/  I.ii,*  ob,  t » <"•(  # IGt'l*  1 ->t  1»  <JSB1»  k.s  ' 1 

i f i ] ! J Mr  (J,  J J 

J / C 0 i / 1 D ; I r t t 1 i 1 1 r P t f < 9 -13*  43 

F - J 

t (J  i 

Lb-  / V *1.4 

BOUNDARY  FACING  l EFT 

it  (.>!.'  ( T »o)  . .!  .1!  . >R.  V ( I M r J)  .’If  . 1 ) GO  TO  1 
r - * + 1 

W b ( i\  ) " i 

GO  T J 4 

BOUNDARY  FACING  FIGHT 

i If  (tt.L  (Irj).  iL.i,  . iv.  !.  )(  I- 1 * o)  .NE.l  ) GO  TC  2 

F - tv  ♦ 1 

1 « j ( n ) r 2 

GO  I l/  4 

BOUNDARY  FACING  UP 

«_  If  ( ,iil  ( I » v ) . jl  . 0 . oK  . r.  J ( I f J- I ) . NF  • l ) GO  TO  3 

F - k + 1 

lVS(i\>  t 0 
GO  I j 4 

BOON  AHY  FACING  COlwN 

•i  i F < :j  m 1 ( I f v I . it  . J . ,ji<  . ‘ lilt  J+l  ) , NE  . 1 } GO  JO  7 

K.  — 1 

IVjlK)  : 4 


4 IMr  «l  f.^lN)  >0  ! 0 t> 

ph in i o 

G F Oki'.mT  ( 1H1 1 / //»  3 AX  « ji'rjuD  CULT.  DIP.  HAS  BCt.'!  F YCEEDE  D ) 

c ah.  t x i r 

t)  IV  (rv)  = i 
J V ; ( f\  ) - J 

K V ii  - K 

7 to  j i nr 
a cO  . i i'  Nf 

C toltlvN  <0o  IDS 

K - J 
LO 

OF-  - . 

O lb  = J 

L>C  1 c.  %J — J F f f l 
uO  To  ( 1 J t i 1 ) » Is 
1J  ltlnl  ( I .Gf-.i  ) uO  1 IP 
t — rv  ♦ I 

IF  (r>.GT,2u  i)  GToP  VUi.  i 
i St  (H)  = I 
uSuOU  - o 
is  - 

GO  I w \s 
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11  lMuWi  ( 1 » J)  ,tv..U)  wO  1 ■ ip 
OSt(M  = w-1 

JB  - j 
00  To  9 

12  CONTINUE 

13  COM  x r ;ur 
KS|j  — K 

C KOV  MOONt'b 

K = u 

00  lo  J — d i i \i  1 

lb  = ? 

14  OS  : 1 

oo  u i = in , r-- 

00  TJ  (lt>»lb)»  j', 

In  1 »-  ( rj.  J ( T * J)  .et.l ) C-.0  1 17 

K = K + 1 

if  (k«ct.2oo)  stop  vb,„;  i 
ISnl (H ) = I 
oSbl  U ) = j 
ob  z c 
00  To  17 

Id  IF  (out  ( I »u)  .k.0. Li)  uo  TO  17 
IStl (K ) = 1-1 
lb  = l 
00  To  14 
1 7 COM  ir.Uf- 
lb  COMl  iUF 
k 5d  1 z y 

Kf  lurt'i 
Lfiu 
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y TT  -J, 
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CO  -.1  iv'  / I iM  T KC  / O'.  6 A » . iFKAL,  bFVBNp,  SFVOKT.  MPR1, 
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b f.OH  = K 

4 OAll  LORD  ( INb.  I .*OKL) 
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NF  KAwF  - huCM  ( U ) 

CALL  LAROlb 
oO  To  2 

9 IF  ( UORptl.L  . *SF\/  1HU*  ) ( TO  10 

SFVBNl  = KNCH(O) 

CALL  VAF^Olb 
00  TO  2 

111  IF  ( I wOHF  »l  .L  • * SF  70H  T * ) iO  TO  11 
jFkokT  - Hnli'-(u) 

CALL  . Ar^OI'-, 

00  To  2 

11  IF  ( 1 *CHr  . »i»PRl  ♦ ) 1.0  TJ  12 


F3PH2 . 

IWHTTF  . 
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I 


\ j 


1 4 


! •> 


1 n 


1 7 


li 


1' 


i ' U 


c 1 


-■V 


C 4 


I f- 1 1 K ( ) 

( Al  w V a U 1 j 

VL  I „ V 

I t ( i >•  • >«l  • . • n 

!>!  » *.  ' I ■ 4c  ( u > 

C uLl.  i A.  U i 'j 
vO  T J • 
it  k i .\ ' t . t 
ct  i,.  _ rl 

C h L l vAm/ 

v O T O L 

ii-  l i .H  Hi  , 
ut  t 4 — :>l 
CALL  VA-OIj 

00  1 o ^ 

Jr  ( J a ( f 
t I-..I  i_v 

1 Ai  L V'  Arn 
oo  1 0 <* 

it  t 1 ...  t- 


' ) A T j li 


. ' 'Hr  5*  ) 1)14 

( ./  ) 


. •-’IVn)  lu  lb 

(u> 


.(  . * i ■ 1 >f  r * ) c TO  16 
AL'"  ( f ) 


. • ,R»  ) vU  T i 17 


I r - .i.l') 


lAlL  . Al-Ui.> 

6 0 Tv  c. 

It-  ( l..ru  uu  To  lft 

L t — t //I'V,  l j) 


CALL  V A 4 iJ  1 ' 

PU  Tv  b 

it  ( i H<  . .L  . tf  v < • ) u To  19 

t V U H “ ' » i V ( U ) 

CALL  LV  Li'. 

lO  To  ?. 

it-  ( I ...  i •(<'■  . ,u  . »r.  STk  * ) T J «)0 

L bl  r\  - ' m.  ( J ) 

l«lL  'A/'Uio 

0 v 1 v ? 

IF  ( .Ur  Kt  .i.t.  • I/MTE*  ) .0  TO  bl 

1 wi<  1 i j z p i A 1 ( C ) 

CAlL  VAADIb 

00  To  ? 

it-  ( U(  HI  .IJL.  '418:.  * ) ,•  T ) <LC 

Olt’K  = -'Nor  (u) 

CALv.  VAf-uib 
Ov  Tv  £ 

it  ( i Al  ;<|  . iL  . •F'LJT  * ) oO  Tv)  t?4 
CALL  l Or o j j ( t L. I A) 
i.t  Lot  - r.HLor  ♦ i 

1 F ( i jPi.OT  • w >»  • 1 ) lAlL  ir  LALl 

C A l i_  1 Z'.'.ii'i  ( t L # i A r . f 1 1. v f » jK  ) 
IF(OK.Lf.i)  Call  T YPl  J . ( N,  AHlUNv 
CAi  l t NCON  ( i.N) 

1A 

t OHM/vT  ( 14,  *&•  ) 

CALL  T rPLlt.i  • UPLo  I Or  a*  ) 

CAi.L  TYPLJM  l';6) 

CALL  TYPLiM*  t,j  vb.iA*) 

CAlL  INTrttT 

IF  (i.vCHf  . -l.  •^L.iTkT  * ) X TO  .b 
HF.  uT  i<  1 “ . ThuF  . 

CALL  I wT  <i_  | 


0 K — 0 a • ) 


2^9 


cb  11-  l 1 .*i  f<!  .!  t; . •Kf.uu*  ) kj  ' T ) £0 
KL^u  '•  . T h OF’  . 

C ALL  1 NT  Kr  T 

20  1 F ( 1 vs C FT  I • ! .t-  • • t Pi_0  1 • ) 0 ) TO  2‘» 

lb/,  — RLU.L  ( Tt.  ) ♦ l*o 

o 0 I u ( 2 7 > i B ) # lb./ 

21  F.PLOF  r .FALbL. 

CALL  1 NT Ki.  I 

2 FI  lF'lOF  = • 1 bl lL  • 

CAlL  JNTRt  T 

29  IF l IftORp.  ,E. »C0n1 • > uO  Tj  30 
LALL  INTKlT 

JO  IF  ( 1 /,  ( rF[.i  • I ,i_  , • V AK  * ) Gu  I 0 3 1 
CALL  VAPOlb 
CO  T v,  2 

31  CALL  I Yl;LI'l(  1 n?U4  1 ) 
oO  T u e. 

SULNUl  T 1 1 A wOKL  ( LNSl » 1 / ORD ) 

Li-ltNMCN  iNbl<22) 

COMMOr  /TtL/  1P01  IT.  ,.c  I»  INb(2L) 
lWCIFu  = * * 

K : u 

00  1 I z \ » Jo  » o 
K = is  + 1 

IF  ( Ii.bl  (K, ) .fih.  .A./,.  IfJbUK  ) .LF-.  *n»  ) GO  TO  1 

1 Pu I NT  z fs  + 1 
Pt  Ti/h; 

1 F Lu(i-1  »6»  lv/OHD)  = f-lI.  ( 0 »6»  1MS11K  ) ) 

I Hu  I , /T  - c 
Kt_  i NK» 

FUNCTION  RnUM(Tu'I) 

UlAL  I.UM2U) 

LOmMJN/TLL/  I HOI  iT  * ,/T  t,  I Mb  ( 2?  ) 

lb//  = 1 
ID  = f, 

II  = 0 

DO  A 1 = JPuiM  * r jCH 
IF  i 1,/  ( T ) .Of  . ’ . « ) Go  1 j i 
I bi,  = 2 
00  1 u A 

111=11+1 

GO  Tu  (’»,_)»  lb// 

2 ID  = ID  ♦ l 

JIJUn(H)  = Il,/S(  l ) - * U * )/?**3U 
4 CONTINUf 

fiUKl  = 0 
DO  b K z 1 » 1 1 
KK  = 1 I + 1 - K 

b uDf’l  = f jDi-I  1 ♦ f/U  -'(KM  4 104  + (K-l  ) 

Htiu//,  z fluaK  fu-iD/flo  T(10*+IL) 

hF.  f UFO  ' 


;?40 


L>Ui  ;• 

Mm  1 M-  , I >(  ,t  k 

■ , I U ■ , !1  RS» EVf  1 

f 

J 9 

» _V0R  * ( ST  Ur 

a ’ , hi  1 Li  , TT I ML  t T V 

HJ 

i.  CVOR  T » Tj  T ■’  # TM  I Sf  »TA, 

I 1 u 

9 1 1 tA  i L> ^ 9 T 7 V 0 1 r 9 T T 7 1 < ND  • T T S T R f N F L ) 

. I ' i S 

1 ON  Sul ( 1 , 1 

( 1 ) » 'h-.M  1)  1 S( 

4 ( 

" 3 » TA  (M  , ?(,)  , 

1 >L  . ( * ) 1 A [ >!  ’ ( 

’ 1 i » t 3 » 1 i ( L u 1 1 * A ) 

f 

ninoi.  jj( ioo) » 

A 1 1 ( 1 9 3 » 1 

( 1 . > l)  ) » i . S ( b U ) t r 

V ( 

h0)»  ITFRSUO) 

VK;/  IlHj.on 

, 9 

L 0 ‘ 1 s. 

/l  'S'/  l AS 

- 1 1 

k t\ 

o 0 ' , ' 

./<  ->-  ./  Mi.  *i  r i 

u >\  i A 

1 / , t tt  ■ ■ * Ik  A » 

Hi . A * r.'l  1 ’ » r.Hk £. / 

i a r a 

i < / t , |0 A * 

t"-' A » u » ?h  I a » 

C-*  1 

T a • ?Hr  a r /"Pa/ 

l.  A I A 

1H4 1 j A 

t jHqUa/ 

041  A 

JO  ul.,  » 

. M'  lA  , . ..  . 0 1 A t 

rj  I 

. 1 a » 

1 . A * 

L t . A < t 1 f.  1 U . A » 

Or 

1 u jn  . a/ 

a ( i 

i t . u i . i ) >u  ru  a 

C «lL 

’ ; 1 1 ( i u • i /• ) 

CALL 

!MTLt.(2) 

c All 

• 1 . T Y t f ( u ) 

CALl 

1 I ' L l r L » , o r I • 3 

CALL 

1 if  1ft  ( } 3 

call 

1 fi(  1 3 » i J o , A) 

CALL 

1 1 1 1 j # it'i  * ) 

CAll 

1 ■'  ( 9 1 # L J t 1 0 # ‘ ) 

LNil 

9(4  il. , .5  3 

CALL 

l r. ( *♦  j»ui»u».  ) 

CALL 

4 ( U / ' 1 i D IJ  Li  9 2.  3 

C 

L N ( Siji  L)C  , j) 

CALL 

I i ( 7 9 C > 'MO*  ' ) 

CALL 

LliCl  • 1 : t • ; 1 J i , u ) 

call 

. 4 'LL  (Poll , rriri  ,1c, 

, T f \ f A M 1 tir,t  T 1C  3/  ) 

call 

l AmLl  ( 6c5 b t .3  J U • 1 .<■ 

cO-.TICiT'r  FtPiCTIor.A) 

c ML  L 

i • ' C c ( 1 ( a 3 1 S kJ  , 1 c 

, i uRAt  ior,^ ) 

C L 

i A • i . i . ( 4 , 0 r b L)  0 » 1 i 

. It-RAT  I or ! A ) 

C A L l 

A ''LL  ( 7m)  r b.,J  > in 

ITLPAr 10SA) 

. • 1 

1 A 1 • ) 

C At_L 

i A 1 ■ i_  u 

, "S'il(l)) 

C Au  L. 

-elf'  , P U - I • 

. • L»0 1(13) 

,0  t 

1 = 1 » u> 

C All 

t Al  Ll  ( ^SS  , J-  1 » 

, "b<o2  (13) 

call 

: f i bL  T ( -i  ) 

o c 3 

JA1 

CO  J 

1 = 1 * A 

11  = 

1 1 1 1 + Uj*I  * ( J — 1 

) *.j0u 

CAlL 

L'J(  1 1 , f • 1 • J » b ) 

call 

i Alii,!  9i’  */  ) 

call 

1 A*  1 l ( I 1-4,  79 » 

:<(!)) 

cor.  r 

i'  tit 

CALL  Kb,  I ( U ) 

CALL  I A>  ti_  (2<i  0 » 1 jL  , .^MLOLaT  lOH  OF  LARGEST  F'PK^R  TFRMSa) 
CALL  <-  F*  bLl  ( 3) 
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Hi  = '0  0 • /FLOAT ( l TLKb ( . ) *5) 

L»0  4 I = 1 . b 
K = 600  4 ( 1-1 ) ♦bu 

CAlL  ( N<  7b J.iWO) 

CALL  LN(  740»K«e.) 

■4  c all  i.  A‘>Fl  ( gbo » * — * i # r, j mu) ) 

LG  b 1 * i * 6 

11  - I 00  + 1 *60 

CALL  l 11(430.  11.0) 

CALL  l.N(4Pu»  II. it) 

b CALL  MJFBlM  3hb.  1 1-4,  ITLi<S(  1 ) *l/b> 

L)0  6 I - 1 » 6 

K = blO.  + FLOAT ( 1-1 )♦ ,2. 8b 
CALL  LU(11u,k,0) 

CALL-  IN(1uU,K,2) 
b L All  1 AilLL  ( 4 0 , K — » » M b u 4 ( I ) ) 
c All  of r sl  t ( o ) 

L uibPl  Ay  !;0UIJL>ARY  PC  If  TS 

LO  9 t-  =1  , f.F  ILL 
Kl  - f LI  ( K ) 

10  7 1 = 1 .h  1 

I ) = xiiP  ( 1 ,K)  - 1 . 

7 Y T L i»U  I ) = YBP(l.K)  - 1. 

LL  = f LMM 

CALL  LUPLT (XTtk, YTLI.Ll) 

11  = > T r M ( 1 ) 

Jl  = TTM4(l) 

CAlL  LN{ 11*10+20 j.Jl *1  i+^0l.)»i) 

LO  b l=?»LL 

11  = >Tf  M(L) 

Jl  = YTfM(L) 

U CALL  I N( Il*l0  + 20d. Ja*!  j + 20G»l) 

9 COuTit  UF 

CAlL  SUI  FLL(FL.IA) 

IAb  = I A 

-.♦♦♦♦***♦**♦**♦#***********  + #***♦*** 

10  lM  = IAS 

Call  F Lf  IUSIFL , l A) 

CALL  LNTYf  L<0> 

CALL  HPL(FL»u»b»IA) 

CAll  LilTYi  L(i» 
ill  = 1 1 n 4 4 * 1 T l i : 

I = yOO.  ♦ ( ALOGlO  (LS(  J TLK  ) > +4, ) *bO. 

ChLl  IU(N1 ,6UU* 3) 

CALL  IN(N1»1»«J) 

1 = oOO.  ♦ Pl*FLOAT ( iTi  Pj( ITtH ) ) 

CALL  LN(i41*3LO*t3(JO#3) 

CALL  l I(Ni+3t 0.1,2) 


1 = u(0.  4 ( ALOG 1 ) ( E.  w ( i T C.R ) ) 4 0 • ) *b0 
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* o 

j ( ..i 

♦ f‘4  9 

J l»Oj 

‘ Ate. 

. - 1 , 1 ♦ 1 1 , | 

f r ) 

. <*  j 

- 1 A 

l »Ll 

l l f-  5». 

1(b) 

1 - 

0 . 4 

( ALO'J 

l J It-.)' 

) ♦ 4 . ) bf> 

0 Al  t 

l Ml  1 

l » i • 2 ) 

cL. 

1 N ( 0 j 

» j . 1 . <■_ ) 

L L 

1 «' ' LL 

(«  b >.  I 

♦ Lo»  o 

' ; <.ma) 

CAl  l 

l A>  tL 

( tl  J > ■ . 1 

♦ * 9 i 

T < • L F i\ A .'LL 

1 Z 

o'  . ♦ 

( .« L 0 lo 

1 1 ( t V 

) 4 2 . ) ♦ br 

c all 

l .(7 

J t 1 , o ) 

C AlL 

N(ObO  » I •?  ) 

0 >>L  L 

L *4 1 : L 

( )b.  I ♦ It. . t. 

- RMOKA) 

C All 

L A • L L 

1 8 7 U » I 

♦ til 

lOLLKAT  Of 

c All 

OFF  u. 

Ho) 

. i*  0 LolAI ION  Of  u.. ‘Oi.i'T  U<MC  • TLWVS 
T ' h K - < • l 

Mi  z 

to  1 A 1-1 , 

.0  11  JZ  1 . i i 

11  1 a - A^AXl  ( r -1  AX  » ( ! A ( I » *J  ) ) ) 

T Ol  — • 2 0 ♦ I A X 

L 0 _ i z 1 » !> 

10  1 tL  Jzl  ,i. 

it  ( A ( T A ( 1 r <j  ) ) . L I . 1 ul  ) ,0  To  12. 

t . - r . , + 1 

1 . (r  < • jl  * 1 u 0 ) Ov.  1 l 1 * 

1 I ( A n ) ~ 1 

k/  O ( 4 i i ) Z o 

12  CO  .1  i\IJF 

Ho'  ITlo  < oO  I AT  T ‘f  I.0CA1I0  . OK  THf  EM«0»  TRVS 


12  K . i z f i - 1 

C A i . i_  l NTYiL  (1  ) 
i 0 i *♦  K z 1 . t*  1 1 
i i ( 1 I ( < > - 1 ) ♦ 1 ‘j  f , u o 
J - ( JJ  ( M -1  ) ♦ 1 0 ♦ «_  ij  0 


14  CAl  L 

i I ML (4  L,  1 , J.  I A) 

Cm.  L 

O F Or  T ( j ) 

CALl 

' A"Ll  ( lbU.fci,  y <i! 

■F 

CALL 

I A'  L>  ( lbU  f Lit)  »K  » » ( 

1 7 

t All 

1 AHLl  ( 1 2>  -i » 'i  i.l  i 1 

, 

' Al  t_ 

L A 1 • Lt  (lb'Jtt'J.flL 

F 

>> 

r 

r 

LAf  L,  (5^0»a0»  1 V(> 

. 

CAlL 

L A ' . L l ( b b U , t. 1 ) , 1 \j  » 

T , 

C .ILL 

l AMLl  (hj'l.'HJ.  IS]' 

i * 

CALL 

L At  Ll  ( 1 1 78 » 2 U , T f- ■'  i ■ 

. 

call 

L A"ll (84  0 , 1 JU . • 1 

r A 

C Ai_L 

1 Dl  It,  ( HbO  , n j . I I v 

1 1 ' 

c «lL 

r UVi»l-  ( Mb') » . . T I v ' 

1 

CALl 

1 UMfJ,  If  -j0»>4.) , ! 1 . 1 

< ) 

• Vf  NO • - . 12»1H&)») 

> liUyor^z  ,1b.  1 H a ) * ) 

■ L ELAPSE  1 TI^Fz  ,F4. 

T TO  CON  j l Mol  z .rb.O 

I ICITY  t OUNiJAMYr  »F4. 
IC  I T Y IIJTrRlOR=  »F4. 

IM  1NTH  I OR:  » F4 , 2 » bit 
l LAl^OUbZ  ,F4,2»bH  S 


P , bH  wINa  ) * ) 
. bH  SEC  A ) » ) 

P , Sm  SF  C a ) • ) 
r»SM  sfca)  • > 
SI  Ca  ) • ) 
f Ca  ) * ) 


l all  mj'm-  u>bU*2u*  n.-.:  ,c  > 
call  F-f-Sul  ( u ) 
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CALL  'jNf)FL.L  (F-L  , 1 A) 

LULL  CFFbt.  I (u) 
ht I Ok f 

♦ *♦*****♦*.♦♦•**♦»♦**♦*♦**♦♦♦*«»»♦ 

•jU-.jKOI  TINL  LM  I A*  I Y*  i‘.  1 > 

COvMOI./KKI/  FL  ( lboo ) * I A 

I F « T l vj  F R f L 

bO  T u ( 3 * C * 1 ) * Tj  1 

1 call  HLANK(l) 

Call  l i Fit.  ( fl  * 1 x * i y * i rt ) 

C Ai_L  mL  Ailtv  ( (> ) 

FLTLK' 

2 CALL  I I'jUFL.  I a*  IY*  Irt) 

F<t  T 'Jl<! 

3 CALL  l f i T V r l ( 1 > 

CALl.  Lint  (FL*  IX*  I Y , 1m) 

F-'L  T UK.  . 


SUbFuuT  I (Jl  L A-3lL  ( IX  * i Y , ICrF) 
COKMOh/FRl/  rL(lbnn).  I A 
I Ml  of  R fl 
CALL  I'- LA  UK  ( 1 ) 

CALL  L I Nl  ( F L * I X * I Y * I n > 

CAuL  LLANK(O) 

CALL  f MAM(FL » ICH* I A ) 
f » L T OK : 


:>LuRbl.T  1 4F  NUMbc-Fv  ( IX  , I V»  jtJV) 

C Ai_L  LHCOuHl'iUF) 

V.RITt.  (23*  1 ) UUM 
1 FORMAT (14,  1HA) 

CALL  I ARtL  ( I A • I Y r I juF  ) 

Rf_  T JKt, 

bUtJKOl  > T I Nl  hU*HK<  IX*  it  , T t.  M ) 
COv.MOr  VFKl/  F L ( 1 KM6 ) * ! A 

INTtbt K Fl 
CALL  l NCOUc ( i uUF  ) 

*K I Tl ( 23  * 1 ) FF  M 

1 F ok v.aT  ( F b . 2 ► IKa) 

CALL  OFF  bL  T ( 3 ) 

CALL  LAhLL  ( Ia*  IY  * 1 1 3U*-  ) 
l\t  TORI 

bUtjROt  TINL  LAtiLl  < I * * i Y , X .FMT  ) 
UlMLNfilON  FMlali  IliUf  (22) 
COi*y,oli/FKi/  FL  ( lb'itj)  * 1 A 
ll.TLbFR  Fl 

CALL  f NCOLL( ItiUF ) 

»M1Tc.(2'**F  T)  a 
CAlL  (LANkU) 

CALL  L It  L(FL*  IX*  1Y,  lA) 

CALL  F'LANr  (O) 

CALL  CHAR  (F  L.  lfiJF  , I A ) 


L0>  oCat-  i1'C«lL 

f U jC  I I O'j  ..CAul  (►  ,i,w) 

L'-lMLNblGf't  F ( t>  1 # J o ) 

h-1  j FHM1.TI0II  KjUUjF  CALCULATES  A MULTIPLICATION 
LaCT'.K,  A POwtJ  i . TO  ut  USLL  In  THE  PM  INTOUT 

blu  - U.O 
l 0 i 1=1  *, 

L-0  1 J=t  «u 

1 t If,  = A^A/  1 ( o Iis  t At  ,bU  < I • J)  ) ) 

IMilv' tATtl  f>C  1 * J? 

SC  al_i_  = 0.0 
KLTUHi 

2 L = if T(AlOU10(oI3) ) 

I.L.»  1 0**L  = jI j 
IF  L 1 j I Of  ( l»L)  «(>T  iu)  t _i  0 T U J 
b C a.  L l.  - lll.O**  l b“D 
00  T u 4 

J SCALl  - ill • 0* ♦ ( 3-u-l j 
4 fctTUKf 
LNN 
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1 F 0 1 ■ l ••  M 1 > » * D A f A 

i'l’t  <■  I ■' > A'H(1U  »•*)»  Ypp(lU0*4) * ULf.(4)f  Mol ’?*>)*  r*-  D(t»J  »2o)  » 

i »i  > ( i.  o u » -I ) » i ,(ioo»4),  i !<f  ( l o o » 4 ) » jnr(iou*4), 

i S(1(M*‘,),  Jf ' j ( 1 Uu  » 4 ) * NKiPF  (h)  t UJPS(b).  II  (’(2) 
i dTl of  A >1  .l-  » LA  ir.  * 4 if  I L 
LCH.1L/  L ii.  f I c < 

i *v  I A ) L A i Jh / fall  / 

u'Alu  1 L"'F  1 i_/  4 ‘P,,  A/ 


M A J It  lfi  (Hi  IFF  ) 

1 FOf'-  AT  (13a*  A'4  t I OX  * A . ) 

iFSl  - .21 

lI  ic  1 . - LF‘j1 

L.U  c I - I 4 U 1 
i O »'=1  »3i 
. - O 2.  4 - ] * 4 

2 A ( 1 » vJ ) r j . 


Call  AT  A 1 ( OF  1 Ll  * X FiY  f',ITLfO 
FFi  i\  j 

.}  r Of  A T ( 1 il  1 t 2>uX,  1 1jH,  >•  if.uAKY  FCI  TS//rbX»  '*  ( 1 1 X » 1 MX  , 6X*  1HY  ) ) 
F K ~ 

L 0 4 4 - 1 * I.F  i Lt 

4 KK  r -AXO  (KK  t 1 r j T { nUi . ( • ) ) ) 

DO  *2  T-1  tY  r. 

4 P’K  1 1 4 I 6#  it  (XF.F(i»K)f  YllP(i*K)*  K = 1 r I F I L f.  ) 

. > i t J A T ( 1 1 1 t c 1 4 X t IF*  't  ( F a * F o • 2 i 3X  * F b • 2 ) ) 

i - 0 

I.  - o 

, j / f r l * f,F  I ll 
.(  i z LN  ( K ) 
uO  / I - 1 * l A , 

= xn  (N  * if  ix  ( XtJ1  ( L * K ) ) ) 

7.1-  l •...  XC;  (I..  14  I X ( fhP(l.  tY  ) ) ) 

Ml  - I - 1 

1.1  - ' - 1 

V. 

c CALL.  LATt  i‘  jo  a<kay 

Tf.jl  F '•«  FlolU  Jji'f..;-  .Ihc>  0,4  I OTTOM  AMP  I OP  *0W 

V 

11. 1 - inOO 
111  - U 00 
Ili«i  - 0 
ilc?  = tl 

00  9 )--l*I.Fl_t 
L - ol  NIK) 

LO  9 rK=l»L 
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L = 1 

Lb  DO  b*4  LL=LirLb 
XbtLU  r.  aL(uL  + ,') 

L4  Yb(..u)  r ) ,IlL(J 

GO  To  (4b»b6)»  i Pj 

bb  IhlK.Ll.UL  .Au(J.  <(M  .ul  .FLUAT  (M)  ) GO  70  140 

IPO  = ? 

LI  - l. 

LS  - ! S - . 

GO  7b  bb 

plot  coni  ok  liol.... 

bt,  bl  - ( JJ-l ) * t M l ,,  )/^  4 1 
CALL  I NOLI  (Xr  » Yl  ,wJ) 

C Ai_L  i L(  7 1 (Xt  ( J1  ) f Yr-  < > J ) * b) 

L 0 b^  0 — 2!  • b w 

lh  ( Ib,  .L  T . v ) ijU  10  b / 

ol  Z V b ♦ l-J 

GO  To  b( 

b7  bl  = b 

b A CAlL  PLoTi  (XP(bl ) » 7k(..1)»  2) 
b9  LO.  ,7  if  VJh 

lh  (Lb.Li  .o  ) m 7Di  i 
lb„  - — T b V, 

00  lb  4k 


If. 


I 
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! 0 < ■ ’ V r . *_  v 0 i V l i \>  ( ( 

' M V t.  I-  V ■ l i , t » 'f  Jl  J ) 

0 

1 1 t . i ■ i ji . I < I i (.  'j  ) i ( (Si  * ) 

A f • I " r-v  ;>!  , 

M ■ i_  A- 

■ A \ i,i  — . 

,J  A _ . / V ,_1 Ji\  F ( - - j.  ) 

I «-  = ' . /f-'L'j  A ! ( ; 1 ) 

M - -! 

1 r • -i 


Mi  _ 1.1, 

C t , I i — C,  j(i  ii) 

1 1 i ' i - .>  i , , ( Mi  i ) 

4 4oL  it.-  it  i,Ui  e!,  i#  * ) 

K - , 

i>- 1 1-; ) **  ,.l  r.t!)  jj  f i 

1C  ^ 

00  I „ 2. 

1 K = 1 

t . L 0 t I — p t c^ttL 

- - i\  ♦ 1 

1 1 ( ( ~ !)**,-.  O I • tj  ) ) u i ■ ' 

wi  1 - u + 1C 

1 ji  — I 

uC  1 , * 

A j l — “it. 

IS  - ^ 

4 . C _ : ^,,1 , 

40  ( J t f ,,  ) » 1 j„ 

3 O Z _ 

oC  \.J  7 

t > J - , 1 1 +ol-l. 

7 x I i . (In)  nt.  • u ) <o  !■  , 

— . /* '._  * ( 1 )i ) ( 1 f • 1 M ) — j . ( I r J ~ 1 ) ) 

V - - . j/Mt  * ( S ' ( I + i « J j - if-l  ( 1 — 1 » J ) ) 

I V - J|-  T ( o*o  i-  J * V ) 

AX'  M < »,  A XI  ( i A a A » T v ) 

•<  tO.h' 

' > 1 1 . I 7 » - »tXN  Ao 

->  I On ’ , ! ( 1 1 1 //  Jtx  i '/H.,  f.  ;o-»  Fit..  3) 
V ->0  - -A  . J *l)X  * iCAL-  / .,  • A , 

00  Id  «•* 

!'  - l\  ♦ 1 

If-  ((**1  )**!■. ol.tf)  jO  i 1 , 

>■>1  = . + iC 

1 0 Vi  - 1 

jO  I 11 
IV  wl  r.  - IC 
IS*  _ d 
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1 1 

lC  1 b 

l=ji » r 1 1 » r 

CO  I 0 

(12,13).  i 

I*? 

-J  - L 
GO  To 

L 4 

13 

J r id 

Awl  — L 

14 

1 L (oli 

( I * w ) ,nL . n ) ..j 

T . lb 

Xl  - L 

c ALt  ♦ ( f loa  r c i — i 

) *Ha  ) 

Y 1 - 

CALL*  (1  LOftT  (j)* 

V ) 

CA'.L  ; 

LoT(-Yl,  xl , 3) 

G = . 

/•'£;*  (Si,(  I , jH  ) 

- Ivl(  1*J-1)) 

^ - - . 

b/Ht1  * ( b J(  1 + 1 , J) 

- b'l  ( 1-1,0)  ) 

A.  - 

XI  , V SL  * U 

Y2  - Yi  + 

CALL  i Lr T (-Y«£,  a , 2) 

aJ>  z ,..i  t ,6b W'a*  (CoS.']  *U  - SInHI*V) 
Y3  o Yl  + , d‘)*Y  jC  * ( Cub,  ’ I * V ♦ $li,PI*U) 
CALL  i L O T ( — Y 3 * X 3 * L) 

XL  z.  y 1 + .7‘j+VjtHJ 

Y L ; Y 1 ♦ . 7 o * V o _ * v' 

LA^L  i LOT (-Y4,  m , t , 

lb  CC.  ,T  if  UF 
It,  CO  Fir  IJF 
!<t  TtJHI  * 

Ct.L 
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.1  Fuf  PLOri»F’LLl] 

jU,  jK  Jt  T J i it  f'i.'-1  I 1 ( A » tJ  * r • ) 
OUN.Mor  /PL1  / jCAuc#  f-i- 
COvMOl./r-LTl/  ib*»  1 A i » Ifl 
COMfor  /Ol/  CALC/P*  IY**E 
IMloFR  I i F*t  » CAi_C'-!P 
1 f J Ti_of  k FL  (siUuo) 

IF  ( l ri’L.NL  .CAi-Ci'k)  Go  ’O  1 

X = A*SCAi_t 
Y z t.*r,rALt. 

CAUL  PLOT (XjY.ak) 

00  To  3 

1 IF  (kK.LT.o)  GO  TO  b 
ix  i GC  aLl*  A ♦ UHi. 

IY  - PC»\LL*M  ♦ 400. 
oO  To  ( F*  * . » 4 ) » kK 

2.  IF  ( K a • E.0  • 1 b/  ) oO  1 j i 
CAUL  LiJTYf  l(U) 

CALL  l InE.  ( FL » I X I * IY1. 1 ) 
CALL  LNTYF’E.  ( 3) 
lSw  - <> 

3 Cali  L If'UFL  # IX»  IY#  IA) 

GO  To  b 

4 lb*  ; 3 
1X1  = I> 

I Y 1 i IY 

b KlTuKI 
LNb 
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.1  FO”  LWlT#LWFi_T 

ljUriKUl  1 INi.  Lie'Ll  (x  . Y,r  r T > ) 
IiIVl^SION  * YltLo) 

w 

tPs  - .rBb 

K - 0 

w 

1 K : r\  ♦ 1 

IF  IK**  .GlL.isPl  S)  hlTOH' 

Kl  = H + 1 
K2  z.  k + *i 
d - f ( K ) - Y C > 

L - A(K2)  - X ItO 

C - a ( K ) *Y  (Kid. ) - X(*<l)*Y(K) 
l)L.‘  = Sort  ( a ♦ a ♦ ij*a) 

L-Is  = APS ( ( A* A ( h 1 ) I * Y (K1 ) 
IF (OiS.^T.tPS)  Go  TO  ] 

o tLlMlf-JAft  POINT  hKU,-,  ^ 1 1 iE. 

c 

KHTS  - fPTs  - i 
1.0  i I=H,nP1S 
Ml)  = v(i+I) 

J YU)  = YU-U) 

4 IF ( K 1 .LT.kPTS)  GO  To  i 
kl  r oki  i 


CJ/OLM) 
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LNL 

1 t 0'>  rt‘JuK.1  #w;,.)PLT 

i»Ui>:uL  tii«i  nriDpLi  (XMi'r  y ,p*  i.ll.  mf  ili.) 

Lit  ('  I'-'M  Xl-M  lUU»4)  » TM(’(1UU»4)»  HLri(4),  <TL^(inO)»  YTtM(lOO) 

(.o./t;-  v>r  /ru,'  iv  x it  .*,  r i: 

I’Lul  UPUfjOAk  f l Ji  .F  I v>'  Mb  r I ON  FPO./  OPTICAL  TIjPLAY  INPUT.... 

i 0 j (•  1=1, NF  ILl 
r c’  = ; L<  (|\1  ) 
uO  1 i = 1 » Ki? 

X I , l I ) = aoM  I » K L ) -I  . 

1 = Yut’II.Kl)  -1. 

C'AuL  r L(  Tl  ( X T f MM  ) , rTLtf(l)*  i) 

CA;.u  i NPLT  (Xlf.M#  YTl.  .►  ) 

10  r 1 , a P 

OALL  : L^  TKxrt  Mtts.5)  t Y l f ; ( K«i ) * 2) 

/.  COi.l  i.AJF 
2>  COuTlNUF 

Kt  I UK' 

tNu 


\ 


A F or.'  UUIPOT  »ouipuv 

SUriWuLT  llJt  OUTPUMM,,,,  ,i  , VN » NH<  aME  » DT  , RR  , XMP » YRP , RLM » ' F I LF  * 
1 xDi G » Y Jib » 1HF  • JOF  » \ <PF) 

L U L"bIoM  I Hi  ( 1 GU , 4 ) , JDF  ( 1 GU , 4 ) » NHPF(b) 

Gi  •■j.NblON  i(luo),  Ta ( b 1 » 2b ) » St\(6  l»?b),  AUP  ( f>  1 » 2b  ) » 

1 Vi , ( o 1 , £o  ) , A ( c 1 » «.  • i ) » FL(JOOO) 

L IVtuSlON  ;riP(lU0H)»  JHFMlbUfU),  PLN  ( 4 ) , Y1  (100,4) 

CO.-iNot  /*•  AlN/uWD,  a 

LO.IKOI VOUT/Ib»  IV  I,  II.,  jV,  ITVI  r IVEC  NCLS,  NCL 
CO.lKof/PLf/  bC^Lt » f L 

INTtotP  SliF’LOl,  VNf'LGI,  UPLOT,  v'PLOT , TVPLOT,  f3ND,  FL 


M = r -1 
f.l  = I -1 

L X - 1 . /Fs_OA  T ( I ! V ) 

1 SNPLOT  = IS 
VNPLOT  r iV?J 
UPLOI  = lb 
VPloT  “ lv 
TVPloI  = ITVP 
lVt.Ll  = IVc-C 

If  (SufLof «L^iU)  b(J  7 b h 

2 K =,  j\  ♦ 1 
10  J 1 = 1 ,(>' 

L U .3  J=1  ,i. 

J 7 h ( i , J ) = SM ( I » o ) 

(.All  f Lf'T  (o.o,  0,0»  -M 

C «i_L  SYMHL4(i.U,  2.C*  . .1,  loHSTKFAM  FUNCTION,  0.0,  lc) 

1 S T |<  = 1 
GO  To  24 

4 If  ( VhPLOT  .(_G  .0  ) oG  I vy  t 
h VNPLOT  r 0 
GO  6 1 = 1 • r- 
UO  u J= 1 , i ; 

b 1 a ( , » j)  = vn ( I , J ) 

CAlL  I LOT  ( u . 0 ’ D.O,  ) 

CALL  SYVhL4(  J.0»2.U*O.  I,  U3HVOP7  ICITY  FUNCTION,  0.0,  lb) 
GO  1 o 24 

7 iF(UPlOT.uO.O)  O')  10  lJ 
ci  UPLOT  = 0 

Go  io  1 = 1 , •) 

GO  10  Jr  1 » i i 

If  ( uHi  (I’J).LU.O)  GO  To  ) 

1 A ( I * J ) = 0,0 

GO  To  If 

y T A ( I » J ) = ,b/l  A*(bN(  1 , J4l ) - f,s(I,J-i)> 

1J  lOi.11l.Ur 

CAi.w  ‘-Nr  AkY  ( M » l . # a l.M.ii  F , Juf-  , I j * -.  I * F * T A ) 

CAlL  fLrT(O.i),  0.0,  -)) 

CALL  SYPtlL4<3.0,2..i,u.  1, 10MU  VFLOCITY,  0.0,  10) 

GO  To  24 

11  IF  UPLOT. LJ.o)  GO  TO  1,, 


\ 
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1 1:  Vf’v.jl  - 1/ 

M t o 1 j I r 1 , y, 

L 0 1 L w — 1 » i*  a 

IF  (i  ..l  ( T »sj)  .L‘J.O)  JO  T J in 
1 A ( i » L ) - u • o 
oU  l \j  1 1 

ta(I»v.)  = (!> ,( ,4i  »j)-br:(  i-i  * j) ) 

l b cO:  f ir.ijf 

CALL.  : i ( >t  * Is. » ' i , 1 4 l • i ;*F  » JbF  , rji'.f'F  » T A ) 

CA..L  t'Lnl  (U.U#  O.u#  ) 

CAlI-  SYvul4<3.U.2.i)#u.1»10HV  VELOCITY,  0.0,  in) 

CC  lu  24 

U If  ( I wf’LoT.tG.O)  LO  To  0 

1 7 I </Ll. 0 T ~ (> 

Oo  1 ^ I ~ 1 , 

JO  I-J  jrl»r. 

If  (Ll4(  (!»J)  .L'J.U)  Gl  T ./  If'. 

00  Tu  lc' 

lo  U = -.*>/! A»  (SN(  I » j+i ) - SNlltJ-1)) 

\c-  - . S/LX*  ( 5(4  ( I +1 » <J ) - SN<I-1#J)> 

TA(i.J)  = OCR  T ( T 1 ♦ T 1 + T2*T2 ) 

19  CO.^TlNUf 

C All.  1 Ni 1 Ah  Y ( f*’  » w • ■ 1 » i J 1 , 1 f-.F  » JBF # f J:  lPF  , T A ) 

CALL  I L'  T ( 0 • 0 » 0.0#  -3) 

CALL  brMriL4<3.0#2.U»U. 1, 14HTGTAL  VtLOC I TY  , 0.0#  14) 
ou  lo  24 

2 0 IFiIwtCT.LG.O)  oO  To 

21  IVlCT  = 0 

CA^L  SYvoLc(6.0#  *4.7,  0.1#  15HVEL0CITY  VECTOR,  90.0#  15) 
CAuL  VF.l  V lC  ( f*  » f < » _>N  *bUD  ) 
fL*-,  i SCALt. 

bCALt.  = 7.0/FLOAT  (M-i  ) 

LO  2.3  K=l,NFiLE 
L = jlM(K) 

1.0  22  KK-1  #L 

22  Y1(KK»K)  - -YHP(<h,K)  4-  l.o 
2$  COi.llfJUF 

CALC  iNOHcrm#  XbP#  PLN#  NFlLL) 
bCALL  = TLM 

CALL  f-LoT  ( 9.0  , -2 . 7b » -3) 

M TlJK’J 


24  CALL  i.TMbL.4(3.0»  1.7,  .1,  2HK=»0.0,2) 

CALL  AJOiiLA  (3.40 » 1.7,  .1#  Kf< , 0.0,  0) 

CALL  bYMdl.^lJ.O#  1 . *4 , .1,  lOHFWAMfc  NO.  ,0.0»10) 
F H AML  = NJ-  aAME 

CALL  t!»JMfltrt<4.2»  1.4,  .1,  FHA«!.»  0.0#  0) 

C CaLCI.'LaTl  MAXIMUM  Ai.! i MINIMUM  OF  uKIb  VAl.MFS 

2b  T 1 i T A ( 1 , 1 ) 

T2  = T1 

lO  2o  1 — X • N. 
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f o o~i * i 

IF  (.JhL  ( I*J)  ,J  .1  ) uO  T Jf, 
ri  - a N'  T N x ( T1 . 11(1*0)) 

12  = , MAXl (T2*  T A ( 1 * j ) ) 

26  CO.Ja.'UF 

t 

v.  GLULf-ATE  COUTQwR  LL\,Llc. 

KClI  = MCo 

IFUbTR.to.l)  AiCLl  - I .CLS 
Tb  - ( T?-l 1 ) /FLOa r (NlL 1 + 1 ) 

L'd)  = ri  « ib 

uC  d K=2  f IkCLI 

27  t.(K)  - | U-l)  + FS 

2b  PKihT  2Q*  Tl*  T2*  ( i(t')*  K-I*LCL1) 

29  FORNiaT  ( Ihb*  iuX.  IOHwIm.  GRll  VmLUE=  * F 1 S . 4 . 3X* 

1 llhMAX.  wALUt.-  Fib.  ■»/  lOX  * 14hC0NT0UH  LEVELS/ ( 1 OX  ♦ 4F1S.4)) 


JO  CmlL  )’LOT(l.o»  b.i*  — .' ) 

CALL  IS(  PLT(M*f.*  XOIS,  Y:  I S * T A * liMJ  * A ) 

CALL  ) LOT  ( .u,  <>.,>,  -j>) 

31  CALL  COnTOH(M*  U*  b*  XbIS*  YLiIS.  NCL1*  TA*  IMD»  A) 
CALL  hNDPu  T ( ArlP  * Y.P.  .:Lil*  NFlLE) 

1STR  r n 

CAll  I LOT  (9.0*  -/.!>*  -3) 

IT  ( VUf’LOT.r.G.  1 ) GO  Tu  S 
IF  (UHL0T.(_0. 1 ) o')  TO  h 
IF  ( VT'l  OT.tu.i)  00  TO  12 
IF  (TvPLOl.LO.l)  00  To  17 
lF(lvi  CT.k'u.l)  00  TO  PI 

32  Fl  TOR' 


SUbKoi  TINE  HNDAH  Y ( 1 1* n * 1 * N1  * II3F  » JBF  * NbPF  * T A ) 
blVEllSlOM  1FF(10U»4)»  jRF(10ll»4)»  NbPF(h)*  TA(61*2f>) 
LO  M KK=1*4 

l = w;.pf(kk) 

IF  (L.FO.O)  Go  TO  37 
LO  3o  K=lrL 
I = lbF ( K » KK  ) 

J = Jt>F(K,KK) 

IT  (j.liE.l)  GO  To  32 
T A ( 1 * 1 ) = T A (1*2) 

GO  Tv  36 

33  IF(J.r.L.N)  GO  Tv  34 
1 A ( 1 . h ) = T A ( I * r«  1 ) 

GO  Tu  3f> 

34  IF  ( 1 . ;L • 1 ) GO  10  3b 
T A ( 1 ♦ J)  = T A ( 2 * J ) 

GO  T 0 36 

3b  IF(I.rtE.M)  60  TO  jc 
TA(y*u)  = T A < • J I 

3b  COuTlf.UF 
37  CONTINUE 
KLIUKt 
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I f-01'  lOOPL  T * I bGPL  l 

su;  Kut.  rrriL  i.oplu  . ..  vp.  yf»  rA»f.iur»A) 

l 1 T i'j(.  Wr  1 ON  jUl  ( 1)  1 » <’  i ) » i(ol»  i*.r> ) 

ul/,UMc  ION  T A ( t>  1 . 2 o ) * t(7M*  V ( 7i> ) » T(7b)»  5 ( 7C  ) * 
Ul.-tL.4i  ION  tJA  ( lb)  . .'A  l 7 ) 

FPAl  “>  1 l Hi  it 
LOGICAL  uc.vr.f  FIRM 
. . I;  l v>iou  fl  ( jiioit) 

C0.«MOI  /PL1  / jCALF.  . Fl 
(.0v,40“  /O  1 / CALC  IF.  I Y ‘L 
INTt-uFK  U(iL>.  F L # I i Ft  » CaLCMP 

t'O  1 1 = 1.1-. 

00  i o=l» u 

1 I A ( I » O ) = APb  (TA(I.J)  ) 

M = r ♦ 1 

02  = 1 4 

I'f  - , 80  *Y  P/FLOA  T ( 1-1  ) 
t/X  = f Y 
X r At  4 la 
A No  = 0 • 

lk  - ( . ?*«ng*jx 

t i»C  ML!  = C.U 
lO  2 0=1  *N 
L 0 2 1 = 1 » h 

I F ( bNL  ( I . J ) • t_U  • 2 ) 30  T ' 2 

/bCAut  = AMAAl  (2SLALL.  Af35  ( T A U . J ) ) ) 

COM  P UF 

L'O  JO  1 1 = 1 . M 

1 = ,4  ♦ 1 - II 

a = x - Oa 

Y = -l  Y 

0(1)  = ANo*X 

V ( 1 ) = -.b*X 
UAI1)  = 0 . U 
W A ( 1 ) = O.U 

JO  = I 

l/v  = r,  - i 
00  t)  0=  1 . UJ 
OJ  = jJ  4*  1 
UA(JJ)  r 0.0 
VA(JJ)  = 0.0 

Y = Y ♦ UY 

LUjj)  = Afto*(X  ♦ Y) 

IF  (tlNl  ( I . J)  .lO.U)  00  TO  3 
if  (u.Jf  < I »M  .00. 1 ) GO  T .)  b 
Y(oJ)  = .M(Y-X) 

IKJ.IQ.N)  00  TO  0 

IF  (buM  1 » J4l ) ,LU.2>  aO  TO  h 

UA(OO)  = .0t6*uY*  ( 1 ,-A  ( l *0+1 ) ) 

VA(JJ)  = ,b*OY* ll .-A 1 1 . J4l ) ) 

JO  =00+1 
V ijJ)  = J(OO-l) 

Y ( jo ) = • b*  ( Y-X  ) 4 T«  ( J . J4  1 ) /2SCALE. 

I'A  ( j J ) = uA(uO-l) 


\ 


TFM(75) 


275 


VA(JJ)  : V«(JJ-1) 

GO  To  h 

3 Vtjv>  = ,l.*(Y-a)  + T A ( u J+l  J/ZSCALE 
UA(JJJ  = ANG*uY*  ( i .-A<  I * J+l  > ) 

V A ( J j ) = .b*UY*(l.-AU,  J+l)  > 

IF  (w.FQ.Nl.)  GO  T J 4 
IF (uuClIt J+l) .Ej.G)  cG  TO  6 
IF  l J+v  .GT,,t)  GO  Tv  G 
IF  (uftD(  I » J-+2)  .LT  ,2)  vO  1J  g 

4 JJ  = JJ  ♦ 1 
U(JJ)  = J(JJ-l) 

V(jJ)  = • b* ( Y “X ) 

UA(JJ)  = UA(JJ-i) 

YA(Jj)  = VA(JJ-l) 

GO  Tg  o 

b IF  (dNC  ( I » J+l ) .tG.Z ) T A ( 1 » J + l ) = 0.0 
V(jJ)  = • L + ( Y-X ) + 1a(I»J41)/ZSCALE 
UA(vJ)  = Ai-.G+'J Y * ( 1 • ” +i  ( r . J+ 1 ) ) 

VA(gj)  = . b*jY* ( 1 , -A { I # J+l ) ) 

IF  ( J.LT  .Nr.)  oO  10  h 
JJ  - JJ  ♦ 1 
U(gJ)  = U(jJ-l) 

V(JJ)  = . b* ( Y-X ) 
c COfjT  IMJF 

JOww  - ,Tf<Ut*. 

IF  (II.GT.l)  GO  TO  G 

PLOf  FIrtbT  LltJL  wIThOJT  F)IL)[;EN  LlNt  ALGORITHM... 
CAj_L  + LOl  1 (l>(  1 ) * \H  1)  . 3) 

1)0  7 J=2»JJ 

7 CALL  PLOTl  (l)(  J)- JA(  J)  * V(J)-VA(J)»  2) 

LiO  o K=?»oJ 
0 7 <K>  = — i 0 , 0 
00  TO  27 

H CALL  PLOT 1(0(1 ) # v (1)  . 3) 

DO  lu  K=2 1 3 

10  OAuL  PLOTl  (U(K)-JA(K)  , V(K)-VA(K)»  2) 

S(l)  : (t(l) 

7(1)  = V(l) 

F iHbT  = .F  ALSt- . 

IF  ( V(4)  ,L7  .T(  J)  ) FI,«ST  =.  .ThUL. 

UO  ^1  K=4»jJ 
LO  lc  KK-1 »N1 

IF 4AuS(S(Kk)-U(k) ) .01. EF ) GO  TO  12 
IF (V (K)-T(KK) ) 11.  1 4 » 14 

11  1F(U0»/N)  oO  TO  1G 

CALL  ) LOT 1 (U(K ) , V (K) , 3) 

JOwi.  = .FalSc. 

00  Tj  d\ 

1?  COfiT ir  UF 
FKinr  is*  k 

13  F0HP.A7  (lHUt  30X»  OiiMtLI’.  . . . » 2HL=I3) 

KLTl-Kf 

14  F1HST  = ,1  ALbt. . 

lb  IF  ( .i+O T . LO*N ) Go  To  1-, 

C A l.  L » LOTl  (U(K)-JA(K)  » V(K)-VA(K)»  2) 
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C-0  lu  PI 

I*>  f.  ] ^ (I(rif)  - I (KK-i  ) ) /(S(KM  - b ( KK-1  ) ) 

If  I;,.)'  (11(1-  )-u(K-l)  ) .J.Fk)  Ou  TO  17 
bE  _ i (v  ) 

1>U  f v l ‘ 

1/  M2  r (W(K)  - V (K-l  ) j / (..  (K  ) - U(K-l)) 

bb  ^ (kP*0(K-1)  - :,l*r-  (KK-1 ) ♦ T(KK-l)  - M(V -\)  ) / {'*2  - Ml) 

i •>  ii  - -1*^3  - l<kk-1)>  -*-  tikk-i> 

if- l ,lKK-l)-f  K.bT.Sb  . b'  . bS.  nT  . S ( KK ) +LR ) f,0  TO  20 
If-'  (i  oO  TO  1 * 

CAul  r‘LOT  1 i Sb » 1 I , 3> 

CAuL  I LOT  1 (U(K )-uA ( K > * V ( k ) — V A ( K ) * 2) 

u 0 a I * - • T R OE  « 

‘jO  I u 2 1 

lb  CAuL  rLoTl  (Sb-Un  ( lb  ) * Tr*  2) 

PU  CmlL  PLOTl (U(K)-uA(K) , V ( K ) - V A ( K ) » 3) 

10  Ki.  — . f ALbt-  • 

PI  CO. ,1  i .1  If 

22  b ( 2 ) r i’(  i ) 

I U c.'*  K - 2 . w 

LVO  C,  U O—K.UJ 

If  (lU.or.blM+Li^)  oC  Tu  24 
,3  OOi.I  1NUF 
24  b ( K ♦ 1 ) = OIJ) 

.,10Kt  Tut  ,-UUlVJM  VALUE  Of-  Ttifc.  V ( K ) LIME  IN  ARRAY  T ( K ) 

1(1)  = \/  ( i ) 

lO  c.1  K-2»fjl 
1 1 , ' ( K ) r -1U.0 
uO  d o Kf “K  > JJ 
1 f l U l f K I “b  ( K ) ) cui  <. : ) * Po 
tb  1 c.M  ( k ) - AMAX  1 ( V (kK ) t r { rs- 1 > * TLM  (K  ) ) 

, O OO.'.I  i,‘  Ul 
27  COuTlMUF 

1 0 20  Kz2,ta 
2o  I ( is  ) ~ T £n.  ( 4 ) 

2b  COriTir.UF 

.JO  COMiMIE 

IK  rtf  E.NL.CALC-'P)  CrtL,.  bNLFLL(FL) 

KL  T Unf 

bUoHOUTINL  PL0T1 (A#tl#K) 

Li  I ML  Mb  I ON  FL  (3000) 

COMMON  /Ul/  CALC  IP*  TYPE 
COkKoI  /PLT / bCAU.  # F u 
C0m^u!./PLT2/  lb..*  IXxr  I Y 1 
INTtOfR  TYPE*  CALCmP,  IL 

IF  (T TPE.NL.CALCPiP)  00  10  31 

CALL  I L^T  ( Am>*kK) 

00  lu  3b 

31  IF (KK.LT.U)  00  Tu  ub 
IX  i 70.0*m  ♦ 300.0 
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1 Y = 70 . U*u  ♦ bO’J.U 
00  To  (32»32»34)  * KK 
32  lK(KK.EO.ISW)  oO  IO  33 
CAt.L  LNTYHu(O) 

CALL  LlM£CFL»lXi#iYl*lA) 
CALL  Li4TYF>E<3> 

IS*  - 2 

3 j CALw  1. IML(FL#  Ix»  IY»  I*«) 

00  To  35 
34  IS*  =.  3 
IXI  = I>' 

I Y 1 = IY 

3b  KtTU^f 
LNl 


